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PREFACE 


This work attempts a complete exposition of the modified Hansen’s theory 
developed by Dr. Peter Musen for analysis of the motion of an artificial satellite 
in the earth’s gravitational field. However, any exposition which lays claim 
to being complete is subject to severe criticism, for the sheer mass of details 
that are involved can never be completely covered in a work of practical 
proportions. Nonetheless, it is the attempt of the author to provide a sys- 
tematic presentation which will begin at a relatively fundamental stage of 
celestial mechanics. It is hoped that in this manner, the exposition may be of 
value to those new to the field of orbit computation and to those whose concern 
is primarily machine programming, as well as to those more interested in this 
particular theory of general perturbations. 

Thu- body of this work wa s presented by the author in a series of lectures at 
the Goddard Space Flight Center of the National Aeronautics and Space 
Administration in October, 1960. The questions and discussions which arose 
in the course of this lecture series were of value in determining what were the 
particularly troublesome concepts and techniques in the theory, and an attempt 
is made to deal with them thoroughly in this work. 

The raison d’etre of this exposition is the recently generated high degree 
of interest in artificial satellite orbit computation, and in the Hansen approach 
in particular. The theory has been in use in the computation of satellite 
orbits since Vanguard I (1958 (3) went into orbit in March, 1958, and is the 
basis of orbit predictions of the Goddard Space Flight Center. Despite the 
important role the theory has played to date, its working is not widely under- 
stood, and it is hoped that this exposition will lead to a greater understanding 
and use of the theory. 

The author, at the request of the Data Systems Division of the Goddard 
Space Flight Center, undertook a study of Musen’s development in order that 
an exposition of this type, beginning at a fairly basic level, might be made 
available to the growing number of those involved in satellite orbit computation. 

It is with the deepest gratitude that the author acknowledges the invalu- 
able assistance and guidance offered him by Dr. Peter Musen, who gave freely 
and amiably of his time, in order to make clear the more subtle points of the 
theory. As well, the author wishes to extend his deepest thanks to Mr. David 
Fisher of the Goddard Space Flight Center, at whose suggestion this work was 
undertaken, and whose direction, assistance, and moral support were instru- 
mental in achieving a finished product. In addition, grateful acknowledgment 
is given to Dr. Alan Galbraith and Mr. Lloyd Carpenter for carefully reviewing 
the manuscript and rendering valuable suggestions. 
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SUMMARY 

A comprehensive description of the Hansen theory 
of satellite orbit calculation, as modified by Musen, 
is presented. The equations of the theory are devel- 
oped in sufficient detail to allow the reader to relate 
them to fundamental laws of celestial mechanics. 
The physical and mathematical concepts underlying 
Hansen’s coordinate system and auxiliary ellipse 
are treated. The disturbing potential function and 
its derivatives are developed in the derivation of 
equations for the perturbations in the orbit plane, as 
well as the perturbations of the orbit plane. The 
method is described for determination of the final 
position and velocity vectors of the real satellite, and 
determination of the osculating elements. Finally, 
a brief evaluation of the theory is presented. 

INTRODUCTION 

This exposition presents a complete development 
of all the mathematical relationships used in 
Musen’s theory of the motion of an artificial 
satellite in the gravitational field of the earth, a 
theory which is basically an application of Hansen’s 
lunar theory to an artificial satellite. Since 
Musen’s development is based upon Hansen’s 
classical work, an understanding of the latter is 
most important in comprehending Musen’s work. 
However, Hansen’s theory does not at all lend 
itself to an easy, clear, and simple exposition; 
quite the contrary is true. 

For the past 130 years, Hansen’s theory has led 
to confusion and controversy in the world of 
celestial mechanics, and for the most part it lias 
been avoided. The difficulty in understanding 
Hansen arises, as Ernest W. Brown (Reference 1) 
expressed it, “partly on account of the somewhat 
uncouth form in which it is given in the Fundamenta 
and partly on account of the very unusual way 


in which the perturbations are expressed.” In 
other words, Hansen’s techniques in solving lunar 
perturbations were extremely unorthodox, enough 
so that many of his contemporaries and successors 
violently disagreed with him. Nonetheless, the 
undeniable fact about Hansen’s lunar theory was 
that it worked, and with a high degree of accuracy. 
Here, it is necessary to inspect what are basically 
Hansen’s methods if we are to understand Musen’s 
final result. Though it is impossible in a work of 
this length to cover all the details of Hansen’s 
theory, it is hoped that by dealing with only the 
techniques incorporated in Musen’s development, 
we will have a sufficiently clear and complete 
perspective on the theory as adapted to the 
motion of artificial satellites. 

In general, Hansen’s lunar theory had six 
distinguishing features: 

1. A fictitious, or auxiliary, ellipse is introduced 
and placed in the plane of the instantaneous 
orbit, i.e., the plane containing the instan- 
taneous radius and velocity vectors. This 
fictitious ellipse is of constant shape, and its 
perigee moves in a specific manner. 

2. The angular perturbations in the plane of 
the orbit are added to the mean anomaly of 
the fictitious ellipse. 

3. The radial perturbations are expressed as a 
ratio between the radial distance of the real 
satellite and that of the fictitious satellite. 

4. The longitudes are measured from a “Depar- 
ture Point” in the plane of the orbit. 

5. One function, IT, is found which expresses 
all the perturbations in the orbit plane. 

6. The theory is a general one which handles 
lunar perturbations of all kinds. 

Each of these six features is to be found in the 
Musen development. Only four changes are made 
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by Musen in his development, but each is ingenious 
and very significant. They are as follows: 

1. Musen has used the eccentric anomaly of 
the fictitious satellite as the independent 
variable instead of time. He has developed 
all his Fourier series expressions in terms of it, 
whereas Hansen’s Fourier series were de- 
veloped in terms of mean anomaly. The 
idea of using the eccentric anomaly was 
borrowed from Hansen’s planetary theory. 

2. The method of iteration is used in developing 
the final series forms, replacing Hansen’s 
method of development into Maclaurin 
series. This change was made desirable by 
the existence of fast computing machines 
which handle iterations rapidly. 

3. Parameters designated by the symbol X are 
introduced. These parameters, which deter- 
mine the perturbations of the orbit plane, 
allow the introduction of a rotation matrix, 
an ingenious development leading directly 
and simply to the final position vector. 

4. The rotation matrix is introduced in place of 
Hansen’s development in polar coordinates, 
obviating much of the cumbersome calcula- 
tion required by Hansen. 

So in the final form of Musen’s development, 
the basic idea consists of introduction of a fictitious 
auxiliary satellite which describes an auxiliary 
rotating ellipse of constant shape and moves in this 
ellipse in accordance with Kepler’s laws. The 
position of the real satellite is determined by its 
deviations in time, as well as in space, from the 
position of this auxiliary satellite. The perturba- 
tions in the orbit plane are relatively large com- 
pared to those o/the orbit plane, and are separated 
from the hitter. They are then determined by 
the single function W for which a differential 
equation of the first order is formed. The per- 
turbations of the orbit plane are determined by 
four interdependent parameters, the X parameters. 

It should be noted here that the sixth dis- 
tinguishing feature of Hansen’s method, given 
above, has significance in the artificial satellite 
theory. Hansen’s method allows inclusion of all 
perturbing forces on the moon. Clearly, the 
forces which disturb the motion of artificial 
satellites can be more numerous and more complex 
than those which disturb the motion of the moon. 
Nonetheless, Musen’s modified Hansen theory 


allows easy inclusion of such forces as those due 
to the gravitational attraction of the sun and the 
moon, solar radiation pressure, the ellipticity of 
the earth’s equator, and the motion of the earth’s 
water masses. It is not inconceivable that even 
the force of atmospheric drag can be included in 
the theory. However, this effect, which is the 
most troublesome and difficult force to deal with 
in artificial satellite theories, is not yet sufficiently 
understood to allow its easy inclusion in the 
disturbing function. 

The theory as developed by Musen is of unique 
and primary significance because it is exact for all 
orders of perturbations. This renders it most 
valuable for the accurate determination of long 
period effects, and allows long term predictions. 
However, it should be noted that for low altitude 
satellites (on which the drag effects are con- 
siderably larger than other perturbing effects) 
the theory gives, for practical purposes, orbit 
predictions good for approximately two weeks’ 
time. 

Musen has recently completed a modification of 
his development (Reference 2) which circumvents 
some of its more troublesome aspects. However, 
the basic approach is exactly that which is de- 
scribed here. The major difference is that the 
perturbations are developed by using the true 
longitude rather than the eccentric anomaly as 
the independent variable. The new technique 
avoids the necessity of “starring” and “barring” 
the potential function in the development of the 
basic perturbation function W. It also allows 
polynomial representation in most places where 
the present development uses infinite series. 
Much of the problem of machine truncation error 
is thus eliminated. The new modification is 
also limited to only those eccentricities for which 
Kepler’s equation can be readily solved. Until 
some replacement of Kepler’s equation is found for 
large eccentricities, this limitation will exist. 
The techniques used in the new modification for 
developing the perturbations and finding the 
constants of integration are exactly those described 
here, and the machine program is roughly equiva- 
lent to the existing one. 

The notation in the following exposition is, 
where possible, the same as that used in Brown’s 
discussion of Hansen’s Lunar Theory (Reference 
1), which was based upon Hansen’s notation in the 
Darlegung (Reference 3). 
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SECTION I 


BASIC NOTATION AND FUNDAMENTAL EQUATIONS OF CLASSICAL 

CELESTIAL MECHANICS 


DISCUSSION 

Before an attempt is made to dissect the Musen 
theory, it is necessary to become familiar with the 
basic expressions which occur in it, as well as the 
notation used. This section lists the notation used 
and then gives the derivations of five equations 
common in classical celestial mechanics, which are 
used in the Musen development. Two of these 
equations apply only in an “ideal” coordinate sys- 
tem (defined as a coordinate system which rotates, 
but so that the form of the equations of motion 
of a satellite is invariant). This is a very impor- 
tant consequence, and one upon which Hansen 
relies. Hansen’s particular “ideal” rotating coor- 
dinate system is an orthogonal one in which the 
X and Y axes are allowed no rotation about the Z 
axis (see Figure 1 which is explained in greater 
detail in Section II) . 

The basic ellipse equations found in this section 
are derived in Appendix A, as are the equations 
of motion in polar form. We start our develop- 
ment with the classical two-body problem, and 
since much of Musen’s development is a vectorial 
one, let us begin here by listing all the vectors 
and angles which will occur. We assume an or- 
thogonal, xyz, inertial coordinate system whose 
origin is at the center of the earth, which in the 
first approximation we will assume to be spherical 
and homogeneous. The xy plane is the earth’s 
equatorial plane, and the z axis has the positive 
direction of the axis of rotation. In the absence of 
disturbing forces, a satellite of negligible mass has 
as its equation of motion 



where 

r=xi+2/j + 2k 

is the position vector of the satellite, r is its magni- 
tude, and i, j, k are unit vectors along the x, y, z 
axes respectively. The velocity vector is 



i+ dy. + dz 

^dt^dt 


k, 


and the acceleration vector is 


d 2 x . . d?y . . d 2 z 
dt 2 l+ dP i+ dt 2 


k. 


In this system of equations, the equations of 
motion have been nondimensionalized by choos- 
ing as a unit of length the earth’s mean equatorial 
radius, and as the unit of time the square root of 
the ratio of the radius to the Newtonian accelera- 
tion at a distance of one radius from a point mass 
(having the earth’s mass). (See Section X.) 

The satellite’s orbit is assumed to be an ellipse. 
The following notation will be used, and is identi- 
cal to the notation used throughout Musen’s de- 
velopment: 

Ei=a = semimajor axis of ellipse, 
E 2 —e = eccen tricity, 

E 3 =co = argument of perigee as meas- 
ured from the ascending 
node, 

£ , 4 =0=longitude of the ascending 
node in the equatorial sys- 
tem, 

2? 5 ==i= inclination of the orbit plane 
to the equatorial plane, 

E 6 =( 7 o=mean anomaly at the epoch 
(i.e., at time t=t 0 ), 

n—wa~ 3,2 =me&n motion (in Vanguard 
units, where standard 
M = l), 

g=(/ 0 -(- 7 j.(<-f 0 )— mean anomaly, 

E— eccentric anomaly, 

/=true anomaly, 

P=F x id-P„j-)-F 2 k=unit vector directed from 
origin to perigee, 

K=P I i+I?J+P z k=unit vector normal to orbit 
plane, 

Q=RXP 

r°=unit vector in direction of r, 
n°=unit vector normal to r, 
lying in the orbit plane. 

In the two-bodv problem, the elements 
E t {i= 1, 2, . . ., 6) are the constants of integration, 
and the complete solution is given by the following 
set of classical equations (Reference 4, p. 164): 


E—e sin E=g, 
r cos /= a (cos E—e), 


r sin 


f=a^l—e 2 sin E, 

a( 1—e 2 ) 
1+e cos / 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


r=a{ 1—e cos E) 
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x, y, z ) 3 sets of 
X, Y, Z > orthogonal 
in, n, Z ) axes. 



i = angle between planes 
f = angle OPr 
Cl) = angle mOP 
v = angle Xr 


Figure 1 . — The coordinate systems of Hansen’s theory. 
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Equation 2 is Kepler’s classical equation relating 
the eccentric anomaly to the mean anomaly. 
Equations 3, 4, and 5 are the standard ellipse 
equations and are derived in Appendix A. The 
final equations necessary to the complete solution 
are those of the rotation matrix and the position 
and velocity vectors: 


-p x 

Qx 

1 


"cos 6 

—sin 6 

0" 

Py 

Qy 

Py 

= 

sin 6 

cos 6 

0 

P , 

Qx 

Px. 


_ 0 

0 

1 . 



“1 

0 

0 ~ 


"cos 03 

—sin ai 

0“ 

X 

0 

cos i 

—sin i 


sin u 

COS 03 

0 


_0 

sin i 

cos i_ 


_ 0 

0 

1 _ 


( 6 ) 


r=Pa(cos E— e) + Q®Vl— e 2 sin E, (7) 


r = -^ (Qa Vl — e ' 2 oos E— Pa sin E)~ (8) 

r\a 

The first and third matrices on the right-hand side 
of Equation 6 each describes a rotation about the 
2 axis. The second matrix describes a rotation 
about the x axis. 

In the two-body problem, the components of the 
position and velocity vectors at the initial time 
t—to can be taken as the constants of integration, 
and the elements can be determined from them by 
means of the above set of equations. The two- 
body problem, however, yields only a first ap- 
proximation to the motion of a planet or a satellite. 
The presence of some disturbing force F causes 
deviations from the simple motion of the two- 
bodv problem, and gives rise to variations in the 
elements of the orbit. Thus, we must find a way 
to deal with these variations. The classical con- 
cept of osculating elements was introduced as a 
device to facilitate the handling of this variation. 
The osculating, or instantaneous, elements of the 
orbit are the elements which would be found at any 
given instant if at that instant the planet or 
satellite were assumed to be traveling in a perfect 
ellipse and in a stationary plane. 

In mathematical terms, the osculating ele- 
ments are defined in such a way that if the position 
and velocity vectors of the orbiting bod} 7 are 
given as functions of the six elements and time, 

x=m,E 2 , . . .,E„t)\ . 

r =g(E u E 2 , . . .,E„t)f’ W 

673 - 619 — 64 2 


then the following relations must hold: 

^dEj dr_ 
hi dt dE t u ’ 

and 

xh dE, dr „ 
hi dt dE t 

where 

s dEj d __ o 
hi dt dE f ~dt 


( 10 ) 

( 11 ) 


is called “Brown’s operator.” (See Reference 4, 
pp. 374-375.) This is a result of the Method of 
Variation of Parameters, a mathematical tech- 
nique commonly used in celestial mechanics; it is 
not a result of deductive reasoning, but rather an 
educated guess as to the form of the general solu- 
tion to the problem. (See Reference 5, pp. 
466-473.) Two consequences of defining the 
osculating elements in this fashion are 


dr 

dt 


. dr — r 
=r and — T > 

d t r 


( 12 ) 


relations which appear throughout the develop- 
ment of certain derivatives of the elements. 

Now that we have introduced a disturbing 
force F which produces variations in the elements, 
we must try to form expressions for these varia- 
tions. To do this it is most convenient to work 
with the disturbing potential 0 rather than its 
gradient, the force. We assume that F has the 
form 


F=grad f2= 


d£2 . . di2 . . dfi , 

,+ ap+aI k> 


dx 


(13) 


which is a convenient form for the development 
that follows. 

In Musen’s development, we require expressions 
for the time rate of change of the angle of inclina- 
tion i, t he time rate of change of the right ascension 
of the node 6, and the periodic part of the time 
rate of change of the argument of perigee g. We 
want these expressions in terms of components of 
the force function. In this development, we will 
introduce the XYZ coordinate system where X, 
Y, and Z are orthogonal axes, the A' axis and Y 
axis lying in the orbit plane, and the Z axis being 
normal to it. Our attempt is to obtain the expres- 
sions in terms of the gradients of the potential 
function along the X, Y, or Z axis. By using a 
vectorial development, we can readily find the 
desired form of these expressions. 
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di 

DERIVATION OF — 

dt 

In the two-body problem, in which we let X 
and Y be orthogonal axes in the plane of the orbit, 
and Z the axis normal to the orbit, we have 



since the force is an inverse square force in the 
direction of the unit vector r°, where r°=r/r. 

If we operate on Equation 14 by rX, we get 

rXr=£xr- 

But 

^Xr=0; (15) 

so, integrating rXr=0, we have 


But we see from Equation 16 that R=A(rXr); so 
R X^=^( rXr)X^ 

=i(rX -)X; 


:[(rA)f-(r-f)l] 


However, we can write 

. dr dr „ , dv dr° 
I= dt=dt r+r dtTt’ 

where drjdt and rdvjdt are the components of r 
along the radius vector and perpendicular to it, 
respectively. It is important here to keep in 
mind that the magnitude of r is not equal to that 
of drjdt, but rather 


rXr=c, 


where c is the vector constant of£integration. 
But from the sketch 



we see that |rXr| is twice the area swept out per 
unit time. Also, the direction of (rXr) is perpen- 
dicular to the plane containing r and dr; so if we 
call the unit normal to this plane R, we see that 


c= 


R 

h’ 


=VQ’+-C-D’ 


So, writing r in terms of its components, as above, 
we have 

dr . dv dr» 

rr=r -3 - t n+rr J,,n' 

which gives 

dr 

T ' T ~ r dt 

because 


dr» 

dt 


and 

Therefore, 


=0 jj - is perpendicular to ^ 


r • r» = r. 


,r h ( . dr\ 

RX ? = ?( rr ^W 


(18) 


where 1 jh is twice the area swept out per unit yy c should again avoid confusion by remember- 
time; and from Kepler s law, hig that in our notation, 


h= ; 

Va( 1 — e 2 ) 

so we have 

rXr=|=c. (16) 

Now, again taking Equation 14 and operating 
with RX, we have 

RXr=-RxJ- (17) 



the radius vector is r, the velocity vector is r, 
and the unit vector along r is r°. Clearly, too, 
r is tangent to the path. That is, dr °jdt is per- 
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pendicular to the radius vector, whereas r is not. 
So, we have from Equation 18 


h 


h dr 


, i iv • iv xjvi 

RXi=i rr — 2 Tut 
r r r at 


But we notice that 


h . h dr 
=- r — jT -r- 
r r at 


h 


d_ 

dt 



Therefore, 


hdr_h dr 
r dt r 2 r dt 



dt 


and we can write Equation 17 in the form 


R X''=—h 


dr« 

dt 


(19) 


Integrating Equation 19 and knowing that 
I* (RXf)=RXr— I* RXA 

where R=0 we have 


RXr+Aro + q=0, (20) 

where q is a constant of integration. 

Now, with Equation 20, known as the Laplacian 
integral, we have a closed vector triangle 



Taking the case where r is in the direction of P, the 
perigee (by definition, the point closest to the 
origin), we see the path is perpendicular to the 
line OP at the point P. Since the velocity vector r 
is always tangent to the path, in this case r must 
be perpendicular to OP, and, therefore, perpen- 
dicular to the vector r lying along OP. Then, if 
r is perpendicular to r, the vector (RXr) must lie 
along r, directed inward toward the origin. So, in 
the case w T here r is along OP, we have both the 
(RXr ) and hr 0 vectors lying along the same line. 
Therefore, in order for the Laplacian integral to 
hold true, q must also lie along the line OP and 
we can write it 

q=AeP, 

in all cases, where P is the unit vector along OP. 
From this it can be shown that e is the eccentricity 
of the ellipse. Now, we know that c=l/h (and is 
therefore a function of the elements a and e) ; and 

c=cR=i R and q=AcP. 


In applying Browm’s operator, 

A— v* d 
dt~~ h. dt dE ’ 

we know that 8/dt applied to any element is equiva- 
lent to the total derivative. So from Equation 16 


which gives 


A —A 

dt C dt 


(rXr), 


dc 

dt 


(jt r ) xi+,x 


Sr. 

dt’ 


( 21 ) 


with hr«, (RXr), and the constant vector q. To 
find the direction of q it is convenient to use the 
following sketch: 



but we know from Equations 10 and 11, 


and 

so w'e have 
Now, 


s t=0 


— r=F 

dt 


dc 

dt 


=rXF. 


( 22 ) 


p dc . dU r N/i 7 i 

dt dt dt dt 


Multiplying through by RX and h, we have 
A(RXR^)+/i(RXc^)=ARX(rXF); (23) 
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but c=l/h, and RXR=0, so 
RX^=MtX(rXF) 


=h(R ■ ¥)r—h(R ■ r)F, 
and since R ■ r = 0, 

RX^=/i ( R-F )r . 

However, R - F is simply the projection of the force 
on the Z axis, i.e., the component normal to the 
plane. Therefore, 

rf=4^ 


„ dR , df2 
RX dt~ h dZ r - 


Since R is a unit vector, dRjdt allows only the 
rotation of the osculating plane about the r vector, 
since r is perpendicular to R by definition. The 
instantaneous angular velocity of rotation of the 
osculating plane, about r, we shall designate as 
ip, a vector clearly in the direction of r. Clearly 
then, 




and we can write 


RX(*XR)=A^|r. 


Therefore, 


(R.R)*-(R.*)R=A ~r. 

But, R- R=1 and R-^=0 since \p is in the direction 
of r, so 

*=A|§r. (25) 

This gives \p, the rotation of the orbit plane with 
the plane considered to be a rigid body. How- 
ever, if we have an ellipse in the osculating plane, 
it is allowed another motion, a rotation in the 
plane, if we disregard changes in the shape of the 
ellipse. If we consider the angle tt between the x 
axis and the line from origin to perigee, then the 
rotation of the ellipse about R normal to the plane 
of the ellipse is clearly (dir/dt)R. Therefore, the 
total motion of the ellipse is given by the rotation 


of the plane in space plus the rotation of the ellipse 
in the plane, and is 

, dfi dw w 

'*SZ r+ S E - 

Considering Figure 2, we can attach three unit 
vectors to the osculating plane: m is the direction 
of the node; K is perpendicular to m and to the 
reference plane; R is perpendicular to the osculat- 
ing plane and to m. Using these three unit vectors, 
we can resolve the total rotation of the ellipse into 
three motions: 

1 . the rotation about ni, which is dijdt, 

2. the rotation about K, which is ddjdt, 

3. the rotation about R, which is du/dt. 

The sum of the three rotation vectors must equal 
the total rotational motion of the ellipse: 

, dfi . dir n di . dO T , , dw „ 

*a z r+ s E -s m+ s K+ s E - (26) 

Multiplying Equation 26 through by m-, we get 

, df2 , . , dir , 

^ 5Z m r m ’ R ) 

=f t (m • m)+f t (m ■ K)+- ^ ( m R). (27) 

But m • m=l, m • K=0, and m • R=0; also, m • r 
= |m||r| cos (/+w) = r cos (»— v). (See Figure 2.) 
Therefore, we have a final explicit expression for 
the first derivative of the angle of inclination: 

h r *S? cos (r— <r)=~- (28) 


DERIVATION OF ^ 

at 

Multiplying Equation 26 through by • (mXR) 
we get 

h r • (mXR)+^ H-(mXR) 

=|m.(mXR)+|K.(mXR)+| R-(mXR). 


But (mXR) is clearly perpendicular both to R 
and to m, so both R • (mXR) and m ■ (mXR) will 
simply equal 0, and we have 

A | |r-(mXR)=^K-(mXR). (30) 
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X 


Figure 2. — Rotation of the ellipse. 


Now, K- (mXR) = m • (RXK). And since R and 
K are unit vectors, RXK= — (KXR) = - (sin i)m 
(see Figure 2), so 

' ( mXR )=-^ sin (31) 

To write the triple product r • (mXR) as a scalar 
function of r, v, and c r, we can temporarily insert a 


set of orthogonal axes x', y', z' along which lie the 
unit vectors i', j', and 1' respectively. Placing 
these axes so that the unit vector m falls along x' t 
and R falls along z ' , we can write: 

r =/’ cos (y — cr)i sin (v — o-)j', 

m=i', 

R=P. 
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Using the standard form of solution of a triple 
product, we have 


r • (mXR) = 


r cos (v— a) r sin (r— a) 0 

1 0 O' 

0 0 1 

— r sin (»— a), 


and finally, 


or 


, bQ , . , S1 dd . . 

h ^ i~~ r sm = ~dt Sln 1 


d6 . . , dfl . . . 

-n sm i=hr sin (v—a). 
dt bZ 


(32) 


DERIVATION OF 


We have from Equation 22 


But, we know that in an ideal system 

A R =°, (37) 

because R, a unit vector always normal to the XY 
plane, does not depend upon the elements, and 
Brown’s operator gives the dependence of a func- 
tion upon the osculating elements. We know, also, 
that since the XY plane is always the plane of the 
orbit, containing both the position and velocity 
vectors, there can be no component of the disturb- 
ing force normal to the XY plane ; if there were, the 
orbiting body would move out of the plane. Thus, 
if we write the force 


where 


r=(F)+||k', 


(¥) ~bX l + W J ’ 


But c=cR, so 

ft R+c sf =tXF - < 33 > 

Multiplying Equation 33 through by • R, we get 

~ R • R+c ^-R=(rXF) • R=R. (rXF). 

But R • R=1 anddR/rft • R=0 since R is a unit vec- 
tor and dRjdt must be perpendicular to it. Thus, 

^=R-(rXF). (34) 

Then, since 

1 , dc 1 dh 

c ~h &nd di~~h 2 di’ 

we have 

^=-A 2 R • (rXF). (35) 

Now, if we take the Laplacian Integral 

RX^+A r# +q=0 (Equation 20) 


and apply Brown’s operator to it, we have 

S„. i , . 5 dr . Sh . , <5r° . <Sq . . 

a Rx s +RX * s +r dt +k < 3f » 


with i', j', k' unit vectors along X, Y, and Z, re- 
spectively, we see that F= (F) because bQ/bZ does 
not appear explicitly in an ideal system. So, 
Equation 11 can be written: 


h . ^5_ dr 
dt* dt dt 


—(F). 


(38) 


Further, because h and q are functions of the 
elements alone, we have 


Sq_dq , 8h_dh 
dt~di Rnd dt~dt 

and, as a consequence of Equation 10, 
5ro 


dt 


— 0 . 


(39) 


(40) 


Substitution of Equations 37 through 40 into 
Equation 36 yields: 

RX(F)+r.g+g-0, (41) 


where we have from Equation 35 

^=-k 2 R.(rXF) = -A 2 F-(RXr). (42) 

Now, considering the unit vector n° which lies 
in the orbit plane and is perpendicular to r, we 
know that n°=(RXr°) and that (RXr) = rn° 
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which we put into Equation 42. Equation 41 
thus becomes: 


RX (F) -roAVF • no+^=0. (43) 


If we write (F) in terms of its components in 
polar form, along the radius and normal to the 
radius, where v is the polar angle, we have 


F 


bn , 1 bn 

-V- r<H — no. 

or r bv 


Putting this into the first term of Equation 43 
gives: 

^ (RXro)+i g (RXno)-AV(F.no)ro+^=0, 


(44) 

and since (RXr°)=no and (RXno) = — r», we write 

^»”-(ig+tM F -n'))r.+g-0. (45) 


Y 



where x is the argument of perigee as measured 
from the departure point X a , and i' and j' are 
unit vectors along X and Y, respectively, we see 
that 

P=i' cos x+j r sin x 

and 


q=A«(i' cos x+j' sin x). (47) 


But F= (F) , and thus 


or 


(F • no) = 


50 

br 


„ , 1 50 

r o . n»4 — — no . no 
r bv 


(F • n»)= 


1 dO 
r bv’ 


Using this last expression for q and the final form 
of dq/dt, we could readily derive the classical 
equations for dJdtQie cos x) an d djdtijie sin x)- 
However, these are not used in the Musen devel- 
opment. 

We now have formed three of the fundamental 
equations of celestial mechanics which are used in 
Musen ’s development. They are: 


since no • n«= 1 and r» . no = 0. Therefore, Equation 
45 becomes 


di . dO , . 

dt~ hr bZ G ° S 


dO A 1 dO i , n . da 
— no — — (l+A 2 r)ro= — ^ 
r bv 


br 


dt 


. dO , dO . . . 

sin % -n=hr sm (v—a), 
dtbZ 


y ( 48 > 


or 

^=(r»XR)^+r. (!+»■) (46) 

In the inertial coordinate system, the rotational 
term w X q would have to be added to this expres- 
sion. 

At this point, we should recall that q is a vector 
directed from the origin to the perigee, and is 
given by 

q=AeP, 

where P is the unit vector in the direction of the 
perigee. From the sketch 


^(roXR)f+ro(I^)g.J 

In addition to these, we will require two relation- 
ships which are classical results; however, they 
are valid only in the ideal coordinate system such 
as the one Hansen devised. 

PROOF THAT cos i IN THE IDEAL SYSTEM 

at at 

It is apparent from Figure i that we can divide 
the rotational motion of the orthogonal XYZ 
system into the rotation components along three 
axes: the z axis normal to the reference plane, the 
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Z axis normal to the osculating XY plane, and the 
m axis along the line of intersection of the two 
planes. 

If we examine the motion of the X axis, we see 
it has a rotation component about z which is 
dd/dt. Similarly it has a rotation component about 
Z which is— da/dt, and one about m which is di/dt. 
Therefore, if we let z°, Z°, and m° be unit vectors 
along z, Z, and m respectively, we can write the 
total rotation velocity of the X axis, to, as 


but 


dd da „ , di 

*>=j t Z°+^7 m°; 


dt 


dt 


dd dd . . ... 

Z ° dt = dt ' C0S t_ ^ n ° sin 


(49) 


where n°=Z 0 X/n°, that is, n° is a unit vector in 
the direction of n, which is normal to m and in 
the XY plane. Note that Z is perpendicular to n 
and that Z, z, and n are all in the same plane ; thus 
we can write 


/ dd . da\ n . / dd A di 

" “U cos '~Tt ) Z, +U 8111 ') m+ dt ”*■ 

This gives the rotational velocity of the X axis in 
terms of three orthogonal components. Now, the 
definition of an ideal system is that the X and Y 
axes have no rotation about the Z axis. There- 
fore, we must have 


or 


dd . da 

dt. cosl -Tt =0 


da dd 

di~dt C ° S 1 ' 


(50) 


This conclusion can also be reached geometrically: 


z 



In our ideal system, we allow rotation of the XY 
plane only about the radius vector r, in which 
case the plane, in some time dt, will be displaced 
to the dotted line. The departure point will be 
displaced to X d , and the line between X a and Xa 
will be perpendicular to both the XY plane and 
the displaced XY plane. Thus, a will be increased 
by da, and 9 by dd. And we will have the right 
triangle 



in which da=(dd) cos i. Therefore, since 
as di^O, we see that 



It is clearly evident that in a non-ideal system, 
this condition cannot hold, for in such a system, 
0 could be held constant (see Figure 1) and the X 
and Y axes could be rotated about the Z axis. 
In this case, dd/dt= 0, but da/dt^O, so 

da ,dd 

lt*Ti eoau 

PROOF THAT IN THE IDEAL SYSTEM r Z, cos i 

dZ O'P 

Inasmuch as the disturbing function il to be 
used is composed of the zonal harmonics of the 
earth’s gravitational field (see Equation 77), 12 is 
symmetric with respect to the earth’s axis (Z 
axis) . Thus, 

Q = Q(r, 4>'), 

where <f/ is the geocentric latitude. 



X 
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If q° is the unit vector normal to r and in the 
direction of increasing </>', bil/bZ will be given by 
the projection of 

1 df2 0 

on the Z axis. The other component of gradient 
12 is in the direction of r which is perpendicular to 
the Z axis and therefore makes no contribution 
to bQ/bZ. Now we can write 

dfi 1 512 
Vv=- YT7 C° s /3, 
bZ r o<t> 

where fi is the angle between q° and the Z axis. 

From the spherical triangle determined by Z°, 
z°, and q° 



where Z° and q° are each perpendicular to r and 
the vectors q°, z°, and r are coplanar, the angle 
j3 between q n and Z # is given by 


Hence, 


cos 13- 


cos i 
cos 4>' 


c)ft_ cos i bQ _ 
dZ r cos <t>' ~d<t> n 

or, if we let i^=sin <f>' , then 


f/i^cos <t>'d<t>' 


and 


dft dff 
r rs=^r cos %, 
bZ oi 


(51) 


SECTION II 


HANSEN’S COORDINATE SYSTEM AND THE AUXILIARY ELLIPSE 


GENERAL OUTLINE OF THE PROCEDURE 

In this section, the rotating coordinate system 
and the auxiliary ellipse will be introduced and 
discussed. At this point it seems advisable to 
discuss the entire problem and method of solution 
in order that the entire procedure be put into 
focus. The purpose of a theory such as Musen’s 
is to allow analysis of the effect of forces on an 
artificial satellite and, from them, predict the 
motion and the behavior of the satellite in orbit. 
In this development, only the zonal harmonics 
of the earth’s gravitational potential are taken 
into account, though other forces could be con- 
sidered as well. We want to be able to predict 
the position of a satellite moving in this gravita- 
tional field, once we have established by observa- 
tion its position and motion at some initial time 
f 0 . From the initial observations of the satellite, 
we are able to deduce an approximation to its 
orbit, which will be the auxiliary ellipse. We very 
carefully determine this first approximation so it 
will have a specific motion. The process pf this 
careful determination is essentially one of separat- 
ing the secular motions from the periodic motions, 
both of which are caused by the disturbing poten- 

673 - 619—64 3 


tial of the earth. Out of this separation process 
come two products: (1) a set of equations which 
defines exactly in time the position of the fictitious 
satellite, and (2) a set of equations which gives 
exactly the relationship between the real and 
fictitious satellites in time. Using these two 
results, we are able to determine the position of 
the fictitious satellite at some desired future time, 
and then the position of the real satellite at that 
time. 

THE COORDINATE SYSTEMS 

Hansen’s first step was to introduce a rotating 
coordinate system. He then defined the motion 
of the rotating system in such a way that the 
equations of motion were invariant in it, i.e., he 
made his rotating system ideal. The coordinate 
systems used by Musen are the same as Hansen’s, 
so only one discussion is needed. 

In the Musen development, a right-handed in- 
ertial orthogonal system xyz is introduced, with 
its origin at the center of the earth, its x and y 
axes in the earth’s equatorial plane, and its z axis 
towards the north pole. Then a second orthogonal 
system, the XYZ system, is obtained by rotation 
through three Eulerian angles from the xyz system 
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(see Figure 3). If the XYZ system is originally 
identical to xyz, its first rotation is that of the A” 
and Y axes through an angle do about the z axis. 
Next, the )' and Z axes are rotated through an 
angle io about the X axis. The third rotation is 
that of the A’ and Y axes through an angle — <r 0 
about the Z axis. After these three rotations, the 
point where the A' axis intersects the celestial 
sphere, named by Cayley the departure point, is 
the point from which all angles in the X Y plane 
are measured (Reference 1, p. 60). 



Figure 3. — Geometry of the departure point. 

Now, with the original position of the XYZ sys- 
tem defined by the angles do, io, and <r 0 , we impose 
such conditions as to render it an ideal system. 
The first condition is that the XY plane be always 
the instantaneous plane of the satellite’s orbit; that 
is, the AF plane always contains the instantaneous 
position vector and velocity vector of the satellite. 
The second condition imposed is that after the 
original position of the AFZ system is defined by 
the angles do, io, and o- 0 , the angular velocity of the 
system, considered as a rigid body, have a com- 
ponent of zero along the Z axis. These two con- 
ditions define the rotating coordinate system as 
ideal, and give rise to the two important relations 
developed at the end of Section I. It is helpful 
to prove this result. 

PROOF THAT THE ROTATING SYSTEM IS IDEAL 


-+F. 


Operating with rX, we have 


rXr=— (rXr)+rXF, 


or 


rXr=rXF, 

since (rXr)=0. But 


( 52 ) 


rXf -B/ rXf “S (rXf) ’ 


where in all cases 


rX r = 


(Equation 16) in which l/h is twice the area swept 
out per unit time. So 

w .. (I /R\ „ // /I \ 1 r/R 

rXf dt\h) R dt\h) + hdt’ 

and Equation 52 becomes 

E s(s)+if- rXF - <53 > 

Operating with R • gives 

R - R ®(i)+ E -w(s)= B '< tXF >- < 54 > 

But RR=1 and R-<fR/c/f=0 since R is a unit 
vector, so Equation 54 becomes 


<t ,/r 

dt 


0=R-(rXF). (55) 

Now, operating on Equation 53 with RX gives 

RxR s(l)+ Rx fG)- Rx(txF) - 


yielding 


Rx ^=A R X( rXF) . 


Expanding the triple cross product, we have 

RX ^ =A[r(R ' F)_F(R ' r)] ' 

But R • r=0, so 


The form of the equation of motion which 
includes the disturbing force is 


BX§=Ar< R .F>. 


(56) 
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Operating on Equation 56 with RX and consider- 
ing that 



,/R 

dt’ 


we have 


rfR 

dt 


=A( R-F)(rXR). 


(57) 


But again using Equation 16, we know that 


so we get 


rXr = 


R 

h’ 



+ (R.F)R=-£+F. 


(62) 


We know dR/dt is perpendicular to R, so (Ry. dR/dt) 
is the vector about which we would rotate 
R to get dR/dt. From Equation 56, this vector 
is Ar(R-F), and so we see that all rotation of 
R is about r, a very informative result. We will 
let this rotation vector be to, so 

co=A(R • F)r. 


If we write F in terms of its components along 
the X, Y, and Z axes, 


and let 


dn dfl 

* 5A + dF J + dZ k 

/pj p I dR p 

(t) +dF J ’ 


we can write 


Since the limitations imposed on the rotating 
coordinate system allow it to have only this rota- 
tion around the instantaneous radius vector, to 
represents the total rotation of the rotating 
coordinate system. The differential operator 
which takes into account the rotation to of one 
system with respect to another is 


Therefore, the equations of motion in the rotating 


system are 

• (lr 1 XX 

r-^+coXr, 

and 


.. d . 


,= ii'+“ Xr 

or 



f =S + $Xr+2(»4)+„X( l .Xr)=-I + F . 

(59) 


F = (F) + ^ k '. (63) 

But R(R*F) is just the component of F along the 
Z axis, because R is identical to k', so 


R(RF) = 



(64) 


Upon substitution of Equations 63 and 64 into 
Equation 62, we have 


..jFr dfi 

r “/< 2+ c>Z 


k' 


= -J+(F)+||k'- (65) 

However, the motion of this rotating coordinate 
system is fixed in such a way that the satellite 
always moves in the XY plane, which means that 
there is, in effect, no force component normal to 
the osculating plane in the XY coordinate system. 
This allows us to consider 


When it is taken into account that in our case 
u=h(R ■ F)r, the equations of motion in the rotat- 
ing system become 


r 


dr 

dt’ 


(60) 


and since 



+toX 


f, 


we have, using Equation 60, 


r 


(Pr 

dP 


+h( R- F) (rXr). 


(61) 


and Equation 14 becomes 

f= S = “£ +(F) - (66) 

Thus, we see that the differential equation of 
motion relative to the system rigidly connected to 
the orbit or osculating plane is of exactly the same 
form as in the inertial system. Although this 
result may at first glance seem trivial or obvious, 
it is not to be expected in a rotating coordinate 
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system and must be shown. The importance of 
the result cannot be minimized for it leads to a 
great simplification of the development. 

THE AUXILIARY ELLIPSE 

We now have two coordinate systems defined, 
the inertial system and the XYZ rotating system 
in which the orbit plane is always the XY plane. 
The standard approach in celestial mechanics is 
to introduce, in the plane of the orbit, some first 
approximation to the real orbit. This inter- 
mediary orbit is determined, and then the devia- 
tions of the real orbiting body from it are deter- 
mined. Hansen’s method was to introduce an 
ellipse of constant shape into the osculating plane 
of the real orbit, with a 0 , e 0 , and n 0 =a 0 ~ 312 fixed. 
A fictitious satellite describes this ellipse as it 
moves in accordance with Kepler’s Laws. The 
ellipse is allowed only one motion in the XY 
plane, and that is a rotation about the Z axis, 
directly proportional to the eccentric anomaly of 
the fictitious satellite. Thus, the argument of 
perigee r in the XY plane is given by 

ir=7T 0 YyE, (67) 

where y is a constant called the secular motion of 
the perigee, to be determined in a specific manner. 
The directions and lengths of the radii vectors of 
the real and fictitious satellites are not identical, 
but differ by the order of magnitude of the per- 
turbations. (The constants y in Equation 67 
and v in Equation 69 have nothing to do with 
the inertial coordinates.) 

The position vectors of the real and fictitious 
satellites are related in space and tune. The 
introduction of the time dimension was one of the 
major causes of controversy among Hansen’s col- 
leagues, though it need not be such a great obstacle. 
The first relationship is that the unit vector along 
the radius of the real satellite, denoted by r°, has 
the same direction at time t that the unit radius 
vector of the fictitious satellite, denoted by r°, has 
at the “pseudotime” 2 . Thus, 

r°(t)=r°( 2 ). (68) 

The second relationship defines the ratio of the 
length of the real satellite’s radius vector, at time 
t, to the fictitious satellite’s radius vector at 
pseudotime z as (1+r). In this ratio, v is small 
and can be considered a “lengthening” or “short- 


ening” factor; r is the radius vector of the real 
satellite, and r is the radius vector of the fictitious 
satellite. The ratio can be written: 

r(t) = (l + r)r( 2 ). (69) 

It becomes necessary, therefore, to know the 
relationship between the real time and the pseudo- 
time, or “disturbed time,” as well as the factor v, 
to determine the position of the real satellite once 
the position of the auxiliary satellite is known. 
Let the difference in times be defined as bz, 

bz=z — t, (70) 

where bz denotes the perturbation of time. We 
can write one further relation between the two 
satellites, and that is that the polar angle v of the 
real satellite at time t is equal to the polar angle 
(J+fl'o+l/AE') of the fictitious satellite at time 2 . 
In this expression, / is the true anomaly of the 
fictitious satellite (see Figure 4), A E=E—E 0 , 
where E 0 is the eccentric anomaly of the fictitious 
satellite at the epoch, and tt 0 is the argument of 
perigee of the fictitious satellite at the epoch. 
Again, y is the “secular motion” of the perigee. 
Both polar angles are measured from the X axis: 

v=f+ir 0 +y&E. (71) 

With the true anomaly of the fictitious satellite 
given by / and its eccentric anomaly by E, the 
motion of the fictitious satellite is governed by the 


usual two-body problem equations: 

rcosf=a 0 (cosE—e 0 ), (72) 

7 sin /=u 0 v 1 — Co 2 sin E, (73) 

7=a u (l— c 0 cos E), (74) 

and Kepler’s equation becomes 

E— e 0 sin E=g 0 +n 0 (t— t 0 )+n 0 bz. (75) 


The “area integral” for the fictitious satellite re- 
tains its usual form: 


where 

Ac= , ] - - 

V«o(l — e 0 2 ) 

(see Appendix A). 

From the above set of equations, the motion of 
the fictitious ellipse in the osculating orbit plane 
can be uniquely determined. The essence of the 
problem, then, is to determine v and bz in order 
that the position of the real satellite can be found. 
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X, y, z ' 

X, Y, Z 


2 sets of 
orthogonal 
axes . 


z 



co - angle mOP 
v = angle XT 

7 t = co + a 


Figure 4. — The geometry of the auxiliary ellipse 
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SECTION III 

THE DISTURBING POTENTIAL AND ITS PARTIAL DERIVATIVES 


THE POTENTIAL FUNCTION IN TERMS OF THE 
GEOCENTRIC LATITUDE 

The earth’s disturbing potential function is 
defined as the negative of the difference between 
the earth’s gravitational potential and the poten- 
tial of a perfectly spherical earth of the same mass. 
Musen writes this disturbing potential 12 explicitly 
to the fourth order in zonal harmonics, where the 
first harmonic is obviated because the origin is 
taken at the center of mass. He gives 

12 (l-3*«)+Js(3*-5 >A 3 )+7 (3— 30i/' 2 +35i/' 4 )> 

(77) 

where i p is the sine of the geocentric latitude, and 
(as is shown in Appendix A) is given by 

^=sin i sin (v—a), (78) 

and k 2 , k s , and k i are the geodetic parameters of 
the earth. 

As is generally true in classical celestial me- 
chanics, it is convenient to develop the perturba- 
tion equations in terms of the partial derivatives 
of the disturbing function. In Musen’s develop- 
ment, these partials should have a very particular 
form, that of Fourier series whose terms have 
arguments containing the eccentric anomaly E 
and the argument of perigee (measured from the 
node) co. This form is required by Musen through- 
out the development, and makes much cumber- 
some algebra necessary, but the end result is that 
all perturbation functions can be easily handled. 

In order to get in this desired form, we make 
several transformations. The first step is to define 
co and 6 0 , the angles which designate the original 
position of the XYZ system (see Figure 3), in such 
a way that the following two equations do not 
contain any constant terms: 

2N=<7 0 + 9 0 -ff-d-2aAE, (79) 

2K—{Tq do c-j-O-j-^yAE , (80) 

where A E=E—E 0 . 

We will now show that N and K are periodic 
only. The quantities <r— < t 0 and 0—d o by defini- 
tion contain only periodic and secular terms. 


Secular terms are those which are proportional to 
time. If the secular terms contained in the sum 
of (<r — vo) and (0—6o) are denoted by —2aAE, 
then N contains periodic terms only. This defini- 
tion of 2aAE determines the constant a. Simi- 
larly, by denoting the secular terms contained in 
the difference of (a—a 0 ) and (0— 0 0 ) by 2 ijA E, 
the constant 1 7 is determined and K contains only 
periodic terms. The determinations of a, y, and 
one additional determination of y, the secular 
motion of the perigee, lead to the development of 
x, y, 2 containing periodic terms only. These 
requirements are equivalent to the requirement 
that our expression for the “perturbation of time,” 
n 0 52, contain periodic terms only. 

Rearranging Equations 79 and 80 we can write 

<7=<7o — (« — j?)A E — (N-\-K), (81) 

0=0o-(a + y)AE-(N-K), (82) 

in which the constant, secular, and periodic parts 
of a and 0 are dearly separated in that order. 
Then taking Equations 71 and 81 into considera- 
tion, we have 

— vo) + (y+a — y)AE-\-(N -\-K). (83) 

We next define the mean values of these three ele- 
ments, denoted by (a), (0), and (to), by the follow- 


ing equations: 

(c) = <T 0 —(ot—y)AE, (84) 

(0)=0 o -(a + y)AE, (85) 

and since co=t>— <r— / (see Figure 4), 

(co) = (ir 0 — (To) + (y -\~oi — y) AE. ( 86 ) 

Considering Equations 83 and 86 we can now write 
for Equation 78 

\p=sin i sin [/+ (to) + N-\-K], (87) 

We now introduce four parameters, which we 


can see by inspection contain the periodic parts of 
the three elements, a, 0, and i, and also the con- 
stant part of i. The physics of the problem allows 
no secular motion of the angle of inclination, so it 
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is of no concern to us. The four parameters 
introduced are 

Ai=sin ^ cos N, 


X 2 =sin - sin N, 


X 3 =cos - sin K, 


X 4 =cos -- cos K, 

£ 

and clearly we have the condition that 

x, 2 +a 2 2 +a 3 2 +a 4 2 =i. 

Now if we expand Equation 87, we have 


( 88 ) 


(89) 


i/'=sin i cos (N-\-K) sin [/+(«) ] 

+ sin i sin (N-\-K) cos [/+ (co) ] (90) 


We perform two additional transformations to 
arrive at a final form of i/ in which the arguments 
are written in terms of E and (co). The first of 
these is to introduce two functions l and m which 
have the forms 


l=pr cos L/T(o>)], (93) 

U'O 

m=— sin [/+(co)]. (94) 

a 0 


Considering that 

r=o 0 (l — e 0 cos E), 
r cos/=a 0 (cos E—eo), 
r sin 7=a 0 y'l — g 0 2 sin E, 
we can write 


(Equation 74) 
(Equation 72) 
(Equation 73) 


which, after further expansion and use of the 
relationship 


% % 

sin i= 2 sin - cos 

— £ 


leads to 


1=2 ( 1 H — V' 7 1 — <V) fos [#+(«)] 

+| (l — v 1 — c 0 2 ) cos [E— (co)]— e 0 cos (w) (95) 


\p — 2(X,X 4 — X 2 X 3 ) sin 

+ 2(X2X 4 +X!X 3 ) cos [/+(&>)], (91) 


leaving only/ and (co) in the argument. 

As Musen points out in his original paper 
(Reference 6), the advantage of using the four 
X parameters is one of symmetry, allowing the 
eventual use of the neat rotation matrix. How- 
ever, the parameters chosen give rise to trouble 
if the satellite has an inclination in the region 
around ir/2, that is, it is a polar satellite. For this 
case the three Hansen parameters p, q, and s 
which still contain only / and (co) in their argu- 
ments are useful. They are defined in terms of 
the X parameters (Reference 6) as follows: 



cos N 
coal? 


2 



( , “ 11 0 r 


sin A r 
osl’ 


s= 7 ? =tan K. 

a 4 


(92) 


This exposition will not deal with this special case, 
but the approach followed is exactly parallel to 
that using the four X parameters. 


and 

m=\ (l+Vl — e 0 2 ) sin [E+( co)] 

-i 

— i (l — y'l — e 0 2 ) sin [E— (co)] — <?„ sin (co). (96) 

Substituting Equations 93 and 94 into Equation 
91, we have for i p 

i>= 2 (XiX 4 — X 2 X 3 ) hi — 2 ^ (X 2 X 4 -(-X 1 X 3 )f, 

(97) 

where m and l have been expressed in trigono- 
metric series whose arguments are in terms of E 
and (co). We will next find such a form for a 0 Jr 
and in the course of the development will show 
how the X parameters evolve into Fourier series in 
E and (co). 


DEVELOPMENT OF -- IN SERIES FORM 

r 

Brown and Shook (Reference 7, page 70) 
describe the development of ( a 0 fr) n into a Fourier 
series in E, but for our purposes only the first 
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power of the expression is needed. First we 
define two functions <t> and /3: 

</>=sin _1 ets or e 0 =sin </> (e a is eccentricity), (98) 

and 

• <P 

„ sm o 

/3=tan |= (99) 

cos- 

Hence, 

l~/3 2 A 2(3 /t AA\ 

cos = r+^ and e ° = i+^' ( 10 °) 

Using Equation 100, we know that 

®o 1 l~h^ 2 

r 1 — e 0 cos E l-}-(3 2 — 2(3 cos E 

Writing cos E=(e iE -\-e~ iE )/ 2, where e is the base 
of natural logarithms, and rearranging, we have 

a 0 (l+/3V g 
r (e iE -0) (l—/3e iE )‘ 

This can be expressed as the sum of partial frac- 
tions in the usual manner, and for convenience, 
since /3 lies between zero and one, can be written 

«o 1 + /3 2 f 1 . pe tE 1 

r 1— /3 2 \_l—0e~ iE 'l~0e iE J 

=\±$[ll+Pe-'°+Pe-™+... 

+pe iE (.l+0e iE +0 2 e 2iE + . . .)] 

= -fz|r g +/s cos U+d 2 cos 2E+ . . .]• 

But from Equation 100, 


0= 


«0 


( 101 ) 


so 

go 2 


Vi- 


-e 0 ‘ 


G 


l + Vl-eo 2 
+/3 cos E+i 3 2 cos 2S+ . . (102) 


With the a 0 /r series written in terms of the 
eccentric anomaly E, and with the X parameters 
eventually developed in E and («), it is clear 
that \p is in the final series form desired. If 
placed in the potential function 12, and with l/r 
written as 

1_1 «o 1 
r~a 0 ~ (1+r)’ 


where 1/(1 + r) is a Fourier series in E and («) 
alone, the potential function 12 and its partial 
derivatives could be given as trigonometric series 
in E and (w) alone. 

SEPARATION OF THE TWO ECCENTRIC 
ANOMALIES 

At this point in the development, we come to 
one of the most difficult operations, one which 
usually is a major obstacle to a clear understanding 
of the theory. The fact is that we have to 
distinguish the E entering into the development 
of the perturbations from the “elliptic” E entering 
into equations derived by using E as a geometric 
angle. This distinction is a very subtle and 
confusing one, and must be handled very carefully 
throughout the development and the computa- 
tional process. As Musen points out in his paper 
(Reference 6), E has the usual geometrical mean- 
ing, when describing the motion of the fictitious 
satellite in its ellipse. However, the perturbation 
expressions are developed using the eccentric 
anomaly E as the independent variable replacing 
time. These two types of E must be distinguished 
from each other because the partial derivative 
of the potential function btt/bE is taken with 
respect to the “elliptic” E. 

The reason the separation is made will not be 
found in the physics of the problem. Rather, 
this separation is a mathematical trick to facilitate 
the development of one differential equation 
which requires only one integration to find the 
perturbations in the orbit plane. The expression 
IF containing the perturbations in the orbit plane 
is an explicit expression of the three variable 
elements, h, e, and x, where X is the argument of 
the perigee (sometimes denoted by tt). However, 
the development of the elements in terms of the 
potential function provides only the derivatives 
of these elements. Therefore, rather than per- 
form these integrations to find h, e, and x which 
are then substituted into IF, it is preferable to 
form dWJdE which is written explicitly in terms 
of the three derivatives. This allows us to find 
IF by the single integration of its derivative. In 
order to perform this process, r and /, both func- 
tions of E which appear in the IF function, are 
considered constant. 

The method used to distinguish between the 
two types of eccentric anomaly is to temporarily 
let the “elliptic” E be F, and to consider F 
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constant. Replacing r , f, and E by p, <f>, and F, 
respectively, in Equations 72, 73, and 74 we have 


p cos 4>=a 0 (cos F— e 0 ), 

(103) 

'cfl 

Sd 

o 

ll 

1-0 

* 

(104) 

p=«o(l — c 0 cos F). 

(105) 


Furthermore m and l, in the actual computation 
of i/', will be written in terms of F and will be 
called m* and l*. The “star” operation means 
that E has been replaced by F in the proper places. 

After the integration to get W, which contains 
E and F, the temporary substitution is dropped 
and the F’s are replaced by E’s. Throughout the 
development, this replacement of F by E is done 
by the “bar” operator; thus, f(E)=f(F). The 
procedure was introduced by Hansen for the same 
reason described here; however, his problem was 
to distinguish between the time t entering into 
7 and / through z, and the time in the elements. 
(See Reference 8, p. 304 or Reference 1, p. 169.) 

FORMATION OF THE PARTIAL DERIVATIVES 
OF THE POTENTIAL FUNCTION 

So, in order to find the partial derivative of the 
potential function which will enter the_develop- 
ment of the perturbation expression dW/dE, we 
take * ls given in Equation 97. In the expression 
used, called ip*, we will have a 0 /~, m, and l re- 
placed by ao/p, m*, and l*, respectively. It is 
important to note that the A parameters do not 
contain F, but are always expressed in terms of E. 
So we could write 


0 *= 




fflo 3 (1 + 

+— — - — ( ®°Y ( 3 d*— 5 d* 3 ) 

+ ^(TY7?(fy (3 - 30 *“ +35 *' ,): (107) 

where *j/* is given in Equation 106 with 
®o 2 


p ^l — eo 2 


(h+P cos F+P 2 cos 2 F 


+/3 3 cos 3 F+ . . .), (108) 


m*=Kl+v / l — Co 2 ) sin [F+(o>)] 


— 1(1 — Vl — Co 2 ) sin [F— (to)] — c 0 sin (o>), 

(109) 

(*=-§(l + -V 1 — Co 2 ) cos [F+(u>)] 

+ 1(1 — Vi — Co 2 ) COS [F—(w)\ — e o cos (o>). 

( 110 ) 


The (1 + r) factor in 0*, as is true of the X 
parameters, is always a series in E and («) ; no 
replacement of E by F is ever made in the series 
for v. The form of this series for v will be de- 
veloped in later sections. In the development of 
dW/dE, the partial derivatives rfdO/d/-) and 
50 IdE are needed. Later on, the derivatives of 
the X functions will contain explicitly the partial 
derivatives dO/dF The partial derivative dO/dl? 
is obtained by the formal differentiation of 0* 
with respect to F and application of the “bar” 
operator 


5 0 dO* 

d E~ d F' 


(111) 


^*=2 (X,X 4 -X 2 X 3 )w*+ 2 


(A^+AAsK*, 


(106) 


where do / is a series in F, and m* and l* are 
series in F and (u). After i p* is placed in the 
potential function, one further replacement is 
necessary to get the final form 0* which will be 
used to find the derivatives. Since the l/r* 
factors in the potential function are found from 
the form 



cio/r is replaced by a 0 /p, and the final form of the 
potential function is: 

4 


The partial r(dO/dr) is obtained very simply from 
the differentiation of Equation 77, and use of 


(~) n= (rr) Y-)Y-Y- (n2) 

W \l + v/ \a 0 / \rj 


This gives 


do :+ 2 / l 

dr 


-$‘(ruJ(v) i<3 - 3V+35 c> 

which is easily seen to be 


(113) 


G73-011)— C4 
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r ^f = -3Q 2 -4fi3-5fi4, (114) 

where 0 2 , 0 3 , and Q 4 are the 1st, 2nd, and 3rd 
terms of Equation 77, respectively. This is 
identical to 


r 


m 

dr 


— 3Q 2 *~ 412£— 5Qt, 


(115) 


where Q*, ft*, and ft* are the first three terms of 
Equation 107, respectively. 

The partial derivative dft/d \p is obtained again 
by the formal differentiation of ft* with respect 
to t/'*, and is then “barred”: 


where 


dO dft* 
d\f/~d\p*’ 


(116) 


do* 

d\p* 


+£(iuJ(!)‘ (3 - 15 *” ) 

+— 5 ( 7 - 7 — Y ( =?Y( — 60^*+ 140i/'* s ). (117) 
a 0 \l + v/ \P / 


At this point, it should be noted that careful 
development of the terms in 0* can result in the 
following Fourier series form for 0*: 


0* — cos [iE-\-2j(os)-{-kF] 

i j k 

+SSSS., 3 , * sin [iE+(2j+ 1)( «) ■ +kF\. 

i j K 

(118) 

The fact that the coefficients of («) are in this 
form is purely a result of the form in which («) 
enters into m*, l *, l/(l + r), and the X’s. It should 
be noted that even though \p and 0 are starred, 
they contain E. This is again a result of the fact 
that not all the E’s are replaced by F’s ; the X’s 
and the 1/(1 + r) term retain their form in E. 
The star really indicates only that F’s are present. 

After 0* is differentiated with respect to F, and 
then operated on by the bar operator to get dfil/dE, 
the form is 

COS [iE+(2j+l)(»)] 

sin [iE+2j(u)]. (119) 

* j 

The partials dQ/b\p and r(dO/dr) will have nearly 
the same form as dll/d E, but the indices of (o>) 
will differ. (See Section VII.) 

Now that we have determined the expressions 
for the partial derivatives of the potential function, 
we are able to begin our development of the IF 
function, the basic perturbation function which 
itself contains all the perturbations of the satellite 
motion in the orbit plane. 


SECTION IV 

EQUATIONS FOR THE PERTURBATIONS IN THE ORBIT PLANE 


Hansen, and correspondingly Musen, have di- 
vided the perturbations of the orbiting body’s 
motion into the perturbations in the orbit plane, 
and those of the orbit plane. In the orbit plane, 
the deviations from the two-body path are con- 
tained in v, the shortening or lengthening factor, 
and n 0 dz, the perturbation of the mean anomaly. 
The factor v can be thought of as containing the 
variation of the elements a and e, whereas the 
angle n 0 Ss contains the variations of the argument 
of perigee (the periodic variations) and the mean 
anomaly of the real satellite. The determination 
of both v and n 0 8z is done by finding one function, 
the IF function, which includes the perturbations 
of all four elements. 

In Musen’s development, the object is to express 
the quantity for dWJdE as a Fourier series in E, 


(a>), and F. Getting an expression of t his form 
involves a considerable amount of manipulation, 
transformation, and algebra; and as a result it is 
virtually impossible to relate the developed equa- 
tions to the physics of the problem. In addition, 
the variable elements themselves are lost in the 
transformations, and the final form might appear 
meaningless. It should be kept in mind that after 
IF is originally defined, all efforts are made to find 
its derivative in a Fourier series with arguments 
containing integral multiples of E, ( 0 ), and F. 

In the process for deriving dW/dE, we will 
obtain an expression for dn 0 bzjdE, which we will 
integrate to get a final series form for n 0 8z. Also, 
in the course of the development, we will produce 
a series form for IF which isolates the terms con- 
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taining sin F and cos F. This series will be used 
in the iteration processes to determine v. 

One more remark would be helpful here. In 
this part of the exposition, we are only developing 
expressions, series expressions, for dW/dE, W, and 
dn 0 8z/dE. The methods by which these series are 
used to calculate perturbations are discussed in 
later sections. The purpose of this section is only 
to show the origin of the equations. It will be 
noticed that the final explicit expressions are in a 
form suitable for iteration ; that is, the W function 
will have W in it. It is to be understood that 
W n+1 —f(W n ), where n is the number of the 
iteration. 

Our first step is to define the II 7 function. We 
do this by setting up an expression for the “per- 
turbation of time” n 0 8z. We take Equation 71, 
which gives the polar angle of the real satellite at 
time t equal to that of the fictitious satellite at 
pseudotime z, that is, 


v=f+ir 0 +yAE, 

where AE=E—E 0 and y is the secular motion of 
the perigee of the auxiliary ellipse in the XY plane. 
Differentiating, we have 


dE 

dt dz dt ^ dt 


( 120 ) 


But we know from our “area integrals” (see Ap- 
pendix A) that 


dv _ 1 
dt hr 2 

and 

dj = J_. 
dz h 0 r 2 

Substituting into Equation 120 we have 

1 _dz 1 dE 
hr 2 dt hoT’^y dt 
or 


dz 

dt 


ho 

~~h 




But we know that 


by definition, and 


ho 


r_ 1 

1 

v«o(l — c 0 2 ) 


( 121 ) 


(see Appendix A) ; and, from Kepler’s law, 
no=a 0 ~ 3/2 , 

which gives us 


^-=cvWl— Co 2 . 


Substituting these into Equation 121, we have 


dz = (ho \ V__ dE 

dt \hj (IT v ) 2 TioVl — «o 2 V'o/ dt 


(122) 


or, since hz—z—t and 


we have 


d&z_dz 
dt dt ’ 


dSz_ i | (h Q \ 1 y (r Y dE 

dt '\h ) (1 + r) 2 n 0 Vl — e 0 2 \®o/ dt 


(123) 


We see that Equation 123 gives the time deriva- 
tive of the “perturbation of time” and is expressed 
in terms of two other perturbations h 0 /h and 
1/(1 +»»), as well as the secular motion y of the 
perturbed perigee. However, we wish to separate 
the higher and lower order perturbations, as a first 
step toward developing the expression in a form 
suitable for iteration. We can write 


1 (l + 2r+r 2 ) 2+2v r 2 

(1 + r) 2 ~ (1+r) 2 + (l + r) 2 + (l + r) 2 

2 v 2 

= -1+ I+r + (THY’ 

which is put into Equation 123 and gives 

dhz = h 0 fh n \ 2 /hp\ v 2 

dt h^\h) (l+rYW (1 + r) 2 

7t=(-Y§- (124) 

'ft 0 yl — go 2 'Vo/ d t 

Now, if we collect the first three terms in a single 
function which we call W, we have 

: ‘-Hr) it,' (125) 


where wc can sec that IT is of the order of the 
perturbations, and 


— =W4-f fl ' n \ y2 V ( ilY — . 

dt \h )( 1+r) 2 w, n 'l —e, 2 \«o/ dt 


(126) 
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But now, since we know that the equation of the 
orbit is 

_ a(l — e 2 ) 


1+e cos / 


(See Appendix A), 


but since 


we have 


hn 2 = 


a o0 Co 2 ) 


where in Musen’s method r=(l + i>)r, we can pro- 
ceed as follows: 

We have shown previously (page 7) that 
q=AcP, where P is a unit vector directed from 
the origin to the perigee, and we know r is the 
vector directed from the origin to the real satellite. 
Therefore, the angle between r and P is the true 
anomaly /. And, by the definition of the dot 
product, we see 


w=- 1— : 


ho 


2 hr 


2r • q 


h Aq&o(1 ?o 2 ) h 0 a Q (l €q 2 ) 


( 129 ) 


r q 


=cos f 


or 


r he 
r • q —her cos /. 
Now from the equation, 


(127) 


Now, as has been previously discussed, for the 
development of the differential equation for W, it 
is preferable to separate the perturbations from 
the elliptic motion. Thus, we will replace the E 
by F, r by p, and / by <f>, so that we have the basic 
ellipse equations written as 

p cos <p=a 0 ( cos F—e 0 ), 
p sin <t>=a 0 \\—e 0 2 sin F, 


=Va( 1-c 2 ), 


and 


p = a 0 (l— c 0 cos F), 


we see that our equation of the orbit can be written p 1 p cos )-j p sin (tto 


hr-\-her cos f=\ > 
h 

and if we divide this through by (l + r), 

Krr,) + iud fc '' cos 

But we have shown in Equation 127 that 
her cos / = r • q, so 


Making the proper replacements in IT, we now 
have our first expression for W. To be consistent, 
we should call this W*, but neither Hansen nor 
Musen has used the IT* notation. Thus, we 
shall use W here, remembering that IT has the 
appropriate E’ s replaced by F’s, and that IT has 
all the F’s replaced by E’s. We have 


W r =-1- 


h 0 


2hp 


2p q 


h 


(luO+IU (t ■> ) -*(LU) 


h hofloiX ej) A 0 tt o(l ^o 2 ) 


(130) 


But- 

and 

and thus 


But since 


1 + r 




1+r 


(r q)=r q, 


hr-\- r ■ q = 


h{ 1+r) 

From Equation 128, we can write 
ff=-l-^+2A 0 / t 7+2/i 0 F - q 


(128) 


q=/(eP = Ae(i' cos x+j' sin x), 

where x is the osculating argument of the perigee 
as measured from the departure point (see sketch 
on page 11), and since 

p • q==Acp[cos (5+XO+2/AE’) cos X 

+sin ( 4>-\-ir-\-yAE ) sin x] 

=hep[ cos (4>+TT 0 +yAE— x)], (131) 

we have the classic equation 

Tl + e cos ($+TTo+yAE— x)"| 

L (i-fo 2 ) J 

(132) 


ti t * h 0 .2h p 

14 l ~h + h B a 0 \ 
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Using Equations 103, 104, and 105, we have for W, 
W= -!^ + (f)(±ZiL 

h \ho/ (1- €q ) 

+(j0 l ilgf ) f( cos F ~ e o) cos (TTo+yAE-x) and 


— v'l — e<? sin /'’sin (Tr 0 +yA.E— x)]. (133) 



(138) 


Expanding Equation 133 gives 


1 L h^KhJ (l-e 0 2 ) 

-(f) i V-,,!) cos ('.+V4E-X)] 

+ [(S)<T^ ! T (e C0S [»"+» 4E -»l-«.)] COS F 

+[-(f ) vT=5? si ” ( '"+» aE - ri ] 8i " F - 


(134) 


which we can write as 


1F=S+T cos E+'h sin F (135) 

if we let E, T, and 'k equal the first, second, and 
third bracketed terms, respectively. We see here 
that E, T, and 'k are functions of E and contain 
no F’s. Now, if we multiply T through by e», 
and add it to E, we will have 


S +<0 T=-l-£+! (136) 

This is a very important relation, to be used in the 
determination of h and ho, and it should again be 
noted that it gives h Q /h and h/h 0 in terms of E and 
(«) alone. 

But for now let us turn our attention to the 
development of the differential equation for W. 
We operate on Equation 130 with Brown’s 
operator 5/dt: 


We also know from the equations of motion in 
polar form (see Appendix A), that the component 
of the acceleration in the orbital plane perpen- 
dicular to the radius vector is 


1 4 ( r i fF \ 

r dt \ dt) 


Written in terms of the disturbing potential (this 
component depends on the disturbing potential 
alone), this gives 


dft _ 1 d / 2 (Iv\ 
rdv rdt\ dt)’ 


but we have defined the “area” as 


so we have 


dt h 


d_ /1\ 
dt \h) dv 


(139) 


From Equation 139, another result appears: 

1 dh dfl 


d_ 

dt 


/1\_ 1 dh dtt 

\h) Id dt dv 


or 


dh 

dt 


—h* 


dn 

dv 


(140) 


Then, we recall that q = heP where P is a unit 
vector, and that 


s4?=<'" xR >“+ t, G+ A! )s 


(141) 


sw 8^ q\ , 2 f -8h 8J 1 

dt 0 dt \h)'h Q a 0 (l — c 0 2 ) \_ P df' dt J 

+ A 0 a 0 ( l-co 2 ) \ji' q+P Tt] (137) 

Since h is a function of the elements alone, and ho 
is a function only of the constant elements do and 
e 0 , we know that 


from Equation 46. Furthermore, we know the 
perturbations in p are present only in the vector’s 
rotation due to the rotation of the auxiliary ellipse 
in the XY plane. This rotation of p around the 
Z axis will be given by 

f-(RX5)v4 < 142 > 
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since y(dE/dt ) is the angular velocity of the ellipse 
in the XY plane. ( llearly, the scalar p is not 
affected by the perturbations; that is, 


Taking Equations 138 through 142 into con- 
sideration, we are able to rewrite Equation 137 as 
follows: 


b\V_ , dfl 2 _/ ,, dfl\ 

dt 0 dr h 0 a 0 (l— e 0 2 ) p \ dr/ 

+ «ha (R Xp>V| 

. 2 ldfl_ 0 

^ha 0 ( 1 — e 0 2 ) r dr p f 

, 2fe 2 £>0 _ „ 

A 0 a 0 ( l — «o 2 ) do p T ’ 

and since h 0 2 a 0 (l— c 0 2 ) = l we can write (noting 
that SW/dt=dW/dt ) 


dW = , d<> T 1 _ 2/i 2 

dt 0 dr L V«o(l — < 


,2 _ p- r 0 2/rp ■ r' 

' r hd 2 a 0 (l—en 2 yh 0 2 a 0 (l 

_i 2 /i q — WvD'i ^ 

/io 2 a 0 (l— e 0 2 ) p dr 

i 2 / . dE 

+WUS (lx '»"'I 


f-r° “I 

-e 0 2 )_ 


or, finally, 


dW , dS2 f2p-r° 

dt 0 do I r 


r2p-r° /2A 2 \ p-r°— p ~ 

L r ^W/a 0 (l-e 0 2 )_ 


+^.(,. X R)“ + »l^,g. (146) 

We wish to simplify this expression further. To 
find djj/dt, note that p=pp°. From Equations 103, 
104, and 105, we obtain by differentiation and 
some algebra 

^=a 0 e 0 sin F, (147) 


</<f>_ao-y 1 — Cp 2 

<ZE p 


From these, and the fact that dp°/d<#>=RXp° 
(by the definition of the derivative of a unit vector 
in a plane), we get 


dp^-dp° d<(> . „ dp 

dF' p d<j> dF +p dF 

=^(RXp) a ^ ll _~ e ° 2 


-\-p 9 a 0 e 0 sin F. (149) 


Now, differentiating Equation 130 


w =- '-r+(l)^ 

with respect to F, we have 


e 0 2 )^ h 0 a 0 (l — e 0 2 ) 


dTT 2h dp . 2 dp 

d F h 0 a 0 (l—e 0 2 ) dF' h a a a {\ — <? 0 2 ) dF 

2ha 0 e 0 sin F 2 

/lotto ( 1— ^o 2 ) + h 0 a 0 (1— c 0 2 ) ^ 

• [(RXp) |° y / r=co 2 +p 0 ttoc o sin e]- (150) 

This last equation for dlT/dE can be rearranged 
and written 

2q • (RXp) TblT 2 ha„e u sin X 

ph 0 (l — Co 2 )” [_dX Mo(l— c 0 2 ) 

_2q • pV^o sin XI I 1 

ho«o( 1 f- o 2 ) J \ 1 — e 0 2 

Multiplied through by p and divided by a 0 , this 
becomes 

2q • (RXp) 

A, 0 tt 0 (l— e 0 2 ) 

rpdW . t ,/2/ i p + 2q-p\'l 1 

L«odX f " SU ‘ ( /(„ao(l-Co 2 )/J Vl — e 0 2 

This, with the aid of Equation 130, can be written 
2q • (RXp) 

/lotto (1— Co 2 ) 

=[-- sin X(V+^+lY| - t =L = - 
L«o dE V h /J yi — e 0 2 

Substituting this into our last equation for dW/dt 
(Equation 146), we have 

dW_ h dll / 2p • r° /2A 2 \ p • r°— p \ 

dt 0 dr V >’ + WJa a (l- e( ?)J 


+ 2/i 0 p- (r°XR) ■yy'H 7== 

^ 1 i 
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Now taking our original definition of W (Equation 
125) we see that 

Ttt | I ^o(l+r)+2A 0 A 0 (l — v) 

+ h{ 1+r) ~h( 1+v) 

and 


Putting this into our last equation for dbz/di 
(Equation 126), we have 


f -"’+(£) (ir+ 1) 

V0 


W 0 V 1 to 


(1 + r) 2 
2 (IE. 
dt 


If we multiply through by n 0 and rearrange, this 
becomes 


dngbz __n a W {1 — r 2 )+» 0 r 2 (fT+l) y /rV did 

dt 1 — u 2 yl — e-o 2 \®o / df’ 

or the generalized Hill formula 


dn 0 8z 

dt 


=n 0 


W+v 2 


y 


1 — V 2 


m «*» 


VI— Co 

Now, differentiating Kepler’s equation 

M=E Co sin E=g () -\- Ti( > (t — 1$) -\-Ti(fiz 

(Equation 75) 

with respect to E we have 


1 — c 0 cos E=n 0 


dt dn„8z 
dE' dE ' 


But since r=u 0 (l— e 0 cos E), w~e have 


w 0 


dt _ r dn 0 8z 
dE a 0 dE 


(153) 


Multiplying Equation 152 by dt/dE and using 
Equation 153 in the result gives 


dn 0 8z , 

(W-\- v 2 \ / r 

dn 0 8z\ y ( r y 

~dE~' 

u-v/U 

<IE ) yi— e 0 2 \«o/' 


which simplifies to 

dn 0 8z / if+r 2 \ iT j r 2 / 7\ 

dE \ + 1-0/ 1 -AaJ 

y / r \ 2 dn 0 Sz fl+W\ 

VhvW dE \1 v 2 ) 


Therefore, 


driifiz W+ v 2 / r j 1 — v 2 y { r \ 

dN ~ 1 +W VW Vl flT/ Vl-c 0 2 V«o/ 


(154) 


This is close to the final form of the derivative 
which we must integrate to find the perturbations 
of time. We will discuss its properties later, but 
for now we must use it in the further development 
of dW/dt. We next substitute it into Equation 
153 to get 


; dt r W-\-v 2 ^r ^4_/l— v 2 \ V ( r \ 2 

/ ' il dE~a 0 ~ I f IT W + \l+W vT^cT 2 W 

_ r A W+r 2 \ /l — r 2 \ y /7V 

~<h \ i+Wv +l li+Wy Vw^Vo 2 W ’ 

or, finally, 


dt r( 1 — r 2 ) 

na dE- ao (l+W) 


0 


f ?/r — \ 

O 0 \l —Co 2 / 


(155) 


This is an important equation which will also be 
used in the development of the derivatives of the 
A parameters. 

Now, as a final step, we must express our 
derivative in terms of E. We know that 


dO_dO dr . dO 5/ 
dE dr dE ftf dE’ 


where, by differentiation of Equations 72, 73, 
and 74, 


d/ = aoVl— to 2 
dE r 


dr 


and £ J j£ =a oto sin E. 


And given v—ira +/, we have 

50 dO 

dr _ d/ 

Also, since r=(l + r)r, 


dO _ dO 
r -<—=/• 
dr dr 


So we can write for Equation 156 


dO . „rdfl . dOaoVl — c 0 2 

^=a 0 c 0 sin E _ — tv > 

dE r dr dv r 


dO 


dr (tn-v 1 — en 2 dN Vl 


dO c 0 — sin E dO 

r-—- (157) 


dr 


Taking Equation 157 and our last equation for 
n 0 [dtjdE), Equation 155, we want to substitute 
into our equation for dW/dt, Equation 151. First 
multiplying through by dtjdE, and using Equation 
155, w 7 e have 
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dW_f h r dV fosinE 1 dO~ir 2p-r° /2A 2 \ o f 0 — p lT/ r \ 1— v 2 A. V I.M\ 

dE - \ "Lao/l — «o ii>E V 1 — e 0 2 ^JL r ^~\h a 2 ) «// — e 0 2 ) wW 1 +H'\ VT— ^o 2 °V J / 


SECTION V 


FINAL EXPANSION OF ^ 

a 

We now turn to expanding Equation 15S in 
terms of the two eccentric anomalies E and F. 

If we let 


s =ES : -Kr +1 ) e “ sinf ’]' 


IN TERMS OF E AND F 


r f —h r i dQ r r /2A 2 \ P' r °~P ~| 

— O 7io0 0 dE l_a (A ] — e 0 2 \A 0 2 / u 0 (l— e 0 2 )J 

p,_har_ A rdO r e 0 sin E ~| 

« 0«0 dr l_Vl-e 0 2 J 

P/ _J AoZ. A ^ l~ 2 P • (r°X R) 2e 0 sin E p • r< >~| 

372 r J’ 


_ A 0 r rdfi ! 
noOo* dr _ 


,,, A 0 r , rdO f r 0 sin E /2A 2 \ p ■ r°— p 

-C — A ^ _ l . 2 ) - /I - 2A 

L A'l — Co / 


we have 


n 0 a 0 dr 


ffo(l — «o 2 ). 


dW_ l r A df2 |“2r p • ro r 

dE~ ° Aoao dE L r « 0 / l-r 0 2 a 0 Vl-e 0 2 

r /2A 2 \ p-ro— p ~| 

ffloVl-eo 2 W/ a 0 (l— e 0 2 )J 

A 0 r A rdO T e 0 sin E 2e 0 sin Ep • ro 
+ « 0 a 0 dr l_- v ] -e 0 2 xl-e 0 2 r 

e 0 sin E^2A 2 ^ p ■ r»— p j_2p • (roXR)” 

Vl— e 0 2 \Ao 2 / «o( 1 — e 0 2 ) r 

, S A 


Now, we will deal with rl' alone. Using the 
known relationships 

r 1 

r (1 + r)’ 

n 0 =a 0 - 3/2 , 

A n ~ — J 

Va 0 (l— e 0 ) 

p • r°=p cos (</>— /) = p(cos </> cos/ 


+sin 4> sin/), 


+ (160) r=a 0 (l— e 0 cos E), 

•\1— e 0 

Ap_ «q 

The process of expanding Equation 160 is a long E, /l— e 0 2 

and troublesome one. It is desirable for us to - /, „ „„„ r.-, 

., p=a 0 (l — e 0 cos d), 

expand it term by term, so we write 


^=4' + 5' + C'+Z?'+E'+E'+-tJ!L= 

dE Vl-e 0 s 

where 


p=do( 1 — e 0 cos E), 
r cos/=u 0 (cos E—e 0 ), 


n 0 a 0 dE |_\ r 


B'=h 0 — A^r — 
«0«0 dE L a oA 


j LV r / ooVi — « 0 J 

dn r ? “] 

dEL Un-Jl Cn 2 -J 


> (161) r sin/=a 0 /i _ eo 2 sin E, 

p cos0=a o (cos F—e 0 ), 
p sin </>=a 0A T — e 0 2 sin 7 '" 

We can write U' in terms of E and F. First we 
have 


A'=k 


T ( ~) ~ r -l ( 163 ) 

6 0 2 ) L\«o/ J 
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The term in brackets in Equation 163 can be 
written 

r/2_FX^°1 /^ ico -_ 7 

LVoo/ a 0 J \a 0 /a 0 

_27 cos J p cos 0 27 sin J p sin 0 

&o a 0 do do 

=2 (cos E— e 0 )( cos E— e 0 ) 


+2-y'l — Co 2 sin Ep 1— e 0 2 sin E 
= cos (E— E)+cos (E+E) 

— 2e 0 cos E— 2e 0 cos E+2c 0 2 
+ (1— c 0 2 )[cos (E— E)— cos (E+E)] 
= (2— c 0 2 ) cos (E— E)+e 0 2 cos (E+E) 
— 2c 0 cos E— 2e 0 cos E+2c 0 2 ; 

and we have A' in terms of E and F: 

1 


A' — A 


dr/,,0 / 1 


a) 


[(2 — c 0 2 ) cos (E-E) 


dE Vl+y/ (1— Co 2 ) 

+e 0 2 cos (E+E) — 2e 0 cos F—2e 0 cos E+2e 0 2 ]. 

(164) 

Next, we take term B' which, using Equation 
162, we can write 


B'= A 


da 0 O 1 


W(I=i?)[~ 1+2<,,cosE 

e<? cos 2E-| 2 ]- (165) 


Now, C’ cun be written 


n>— a da 0 Q W 1 

^ ^ x T7» I o 


dE V (1— c 0 2 ) 

r ^ 


L®oV 1 ^o 2 ®o"V 


2 r p • r°— p~| 

1-+ «o j 


( 166 ) 


Since p • r°=p cos (0—/), we can write the brack- 
eted term of Equation 166 as 


[ H 


= " 2^ 

_Vi — 


r r 


e 0 2 ffio ®o -yl — e 0 2 ®o 

27 7 cos / p 


[cos (0—7) 


-«] 


— Co 2 ) «0 Go 


COS 0 


27 p sin 0 7 sin / 

+ ..~ ■ — — - — ■ - 


«0 «0A 1 — e 0 2 ffloVl— «0 2 

2 7 7 p 


V 1— e 0 2 «o a 0 a 0 Vl— Co 2 

Applying Equation 162, and substituting the 
result in Equation 166, we get for C": 

/v_ . da 0 ^ h 2 1 

^ A N 77 LO 


[2 cos (E— E) — 2 


dE V (1 — Co 2 ) 

— e 0 cos (2E— E) — c 0 cos E+2e 0 cosE]. (167) 

Next, term D' of Equation 161 can be written 

ba 0 U 1 f7 ^ . j+1 


D' = Ar 

D' = Ar 
and so 
D' = Ar 


dr 

(1-Co 2 ) 

da tl 9~ 

1 

br 

(l-€o 2 ) 

ba 0 Q, 

1 


br (1— c 0 2 ) 


[(1— c 0 cos E)e o sinE], 
[c 0 sin E— ~ sin 2 eJ. 


(168) 


The next term E' of Equation 161 becomes, after 
consideration of Equation 162 and using the fact that 
P (r°XR) = P sin (7—0), 

50 P 2 p7 

r-yi— C 0 2 

— V COS (/— 0 ) 


E'—Ar I sin (f-0) 

dr Lry 1— c “ 

2p > t - r 


e 0 sin 


in e] (169) 


r w 1 — Co 2 

For the bracketed part of this equation, we write 

[ ]= 


1 


(1— e 0 2 ) r(l+»0 
1 


[2 p sin (./— 0)V1— Co 2 

— 2p cos (7-0) co sin E] 


a ° [(1 — e 0 2 ) sin E( cos E-c 0 ) 


(1— c 0 2 ) (1+y) 7 

— (1— e 0 2 ) sin E(cos E— e 0 ) 

— (c 0 sin E)( cos E— c 0 )(cos E— c 0 ) 

— e 0 sin 2 E(l— c 0 2 ) sin E] 

which after rearranging becomes 


7=7° {sin E[(l — c 0 2 )(cos E-Co) 


(1 — c 0 2 ) (1+r) 7 

— e 0 (cos E— c 0 )( cos E— c 0 )] 

—sin E[(l— c 0 2 )(cos E— e 0 ) + e 0 (l— e 0 2 ) sin 2 E]}. 
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But sin 2 E= 1 — cos 2 E, and a Q fr = 1/(1 — c 0 cos E), 
so we have 


77 — 2 . Trfl 2[sin E( cos F—e 0 ) 

(1 — e 0 z ) (1+x) 

—sin F( 1 — «o 2 ) cos E) 

=rr J 7T^ 2 R sin {E+F) 

(1 — Co ) 1+x L2 

sin ( E—F)—e 0 sin E 

- (1 ~ 2 e ° 2) sin (E+g)- * 1 "* 0 ^ sin (E-E)l 
Thus, Equation 169 for E' becomes 

E - Ar ^[(ub)(lio( e *' S i" ,F+E) 

— (2— e 0 2 ) sin (F—E)—2e 0 sin E)J- (170) 

We next determine terizz F' in terms of E and F: 
first we rewrite F' as 

F> _ . dfflpO 1 27 c d sin E / p • r°— p \~l 

7 dr A 0 2 (1 — e 0 2 ) L a 0 (1— e 0 2 )\ «o /J 

(171) 


Considering the bracketed term, we see that 
r 1 27 P e n sin E 


a 0 a 0 (1 — e 0 2 ) 


[cos (<t>—f)—l] 


2r p e 0 sin E . - . . - . -j , . 

= 7- ijv (cos <A cos /+sni 4 > sin /— 1) 

a 0 a 0 (1 — e 0 ) 

which, after some algebra, becomes 

f 0 sin (F—2E)—e 0 sin E-f 2r 0 sin /?. 
Putting this into Equation 171, we have 


F 


da 0 Q A 2 1 
_ Sr V0=O 

[e 0 sin (F—2E)—e 0 sin F-\-2e Q sin E}. 


( 172 ) 


We now return to Equation 161 which was 

&= A'+B'+C'+D'+E'+F'+- 1 M = • 

VI— e 0 2 

If we replace A', B', C', D', E', and F' by the 
values given by Equations 164, 165, 167, 168, 170, 
and 172, respectively, we can write our final equa- 
tion for clW/dE as follows: 


dW , T . da 0 £2 , , , , ba a Q, 

dE =NAr ~W +MA 


Sy 


Vl — e- 0 2 

where A and S' are given in Equation 159 as 

s= \jt l?“( w,+ I +1 ) e ° sin F } 


(173) 


1-r 2 


_( 1+ L^L A 

1+W\ OoVl — ?0 2 ' 


and M and N are given by 
h 2 

(1— e ( j i )M=-r- 2 [2 cos (E—F)—2—e 0 cos (2 E—F) 

ho 

— e 0 cos F-{-2e 0 cos E] 

+ t( 2— e o 2 ) cos ( F—E ) 

+ C 0 2 cos (E+F) 

—2e 0 cos F—2e 0 cos E+2r 0 2 ] 

— l+2e 0 cos E— i e 0 2 cos 2E— e -^~ 

and 

(1 — e 0 2 )N=~ [e 0 sin (F—2E) — e 0 sin E+2<? 0 sin E] 

fl {) 

+ 1 ^ t(2— e o 2 ) sin (E—F)—2e 0 sin E 
-f e 0 2 sin (E+E)] + e 0 sin E~*^~ sin 2E. 


SECTION VI 

DERIVATIVES OF A PARAMETERS, EXPRESSIONS FOR f-, AND v 

ft tin 


THE X DERIVATIVES 

In Section V, we arrived at the final form of the 
differential equation for dW/dE. We developed 
it in a form such that all its components are 
trigonometric series with arguments containing 
only E, F, and (to). Now we turn our attention 
to the much simpler development of similar forms 
for the derivatives of the A parameters. Several 
relationships we have found before are required in 


this development. 

We have previously defined 2 N and 2 K as 

2N=(To+6 0 —<t—0—2olAE (Equation 79) 

and 

2K=<j 0 — do— o— |-0+2)}AE, (Equation 80) 
and we have developed the expressions 

sin i ( ^n==h r sin (v—a), (Equation 32) 
(it q/j 
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da dd 

dt~~dt C ° S h 
di r dQ , . 

ir h dz rcos (v ~ a) > 

dft dft 

r ^z _ ^^ cos, • 


(Equation 50) 
(Equation 28) 
(Equation 51) 


or 


Also, we have 

dt r 1— v 2 /, , y r\ . 

n 0 ~TEi= I 1-j — — )■ (Equation 155) 

dE Uol+TEV "\ 1 < «• «o/ 

But in Equation 159 we defined 

1 ~\-W \ -y 1 — e o~ a °) 


so we can write 


d_N_ 

dt 


dt _ r 
dE n 0 a 0 

l da 1 d6 dE 


(174) 


(175) 


2 dt 2 dt dt 

Substituting Equations 32 and 50 into this yields 

dN IdO , 1 1 ,)S1 . , , dE 

— — ^ — cos i—^ — . ft ^75 r sin ( v—a)— a 


dt 

or 


2 dt 


2 sin i dZ 


dN dE 1 ,dSi . , v , . . .. 

^=— a ^ — - ft r sin (»— <r)(cot 1,+csc r) 

dE 1 df2 . , , i . . 

= ~ a li~2 h dZ rSm COt 2 ( 1/6 ) 

We now differentiate our original definition 
from Equation 88, 


and get 


Xi=sin ^ cos N, 


</X, 1 i di A7 . i . A7 dW 

dt = 2 COS 2df COSiV “ Sm 2 SlnAr ifF 

from which, using Equations 28 and 176, we get 

dXt 1 i „ /, d!l , A 
dt = 2 C0S 2 C0S N \ h dZ r c0s ) 


+ 


( sill £) 


dE . , T 

a -n- sin A 
dt 


+^sin 0 ^2 r co * \ s ' n sin -W 


dX, 

dt 


. dE 1 , dQ r i A7 , . 

=aX 2 ft ^ r cos - cos A' cos (i v—a ) 

+cos ^sin (v— a) sin wj> (177) 

since X 2 was originally defined as 


X 2 =sin - sin N. 

Recasting the second term of Equation 177 w T e 
have 

dX, dE 1 . dfi 

dJ- aX2 dt + 2 a ° h{1+v) dZ 

I V ^ 

— cos - cos iV cos (v—a) 

L«o ^ 

v % I 

d — cos - sin (v—a) sin N • (178) 

«o Z J 

The bracketed term of Equation 178 can be 
written 

f ]=— cos \ cos [N—(v—a)]. (179) 

do Z 

Since v— (r=/+(co)+W+E’, from Equations 83 
and 86, the right-hand side of Equation 179 is 


— cos 5 {cos K cos [/+ (to)]— sin K sin [/+(“)]} 
a 0 z 

i‘ % 

=— cos [/+(co)l cos - cos K 
do Z 

— — sin [/+(«)] cos ^ sin K. 

But we know from original definitions in Section 
III that 

— cos [7+(w)]=f, (Equation 93) 

do 


dO 


sin lf+(w)]=m, (Equation 94) 


cos g cos K—\ 4 , 


cos - sin K—\ 3 . 


> (Equation 88) 


so the right-hand side becomes 
{X 4 f— X 3 m}; 
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and we have from Equation 178 that 

* <a( 1+0 ||(X 4 Z-X 3 m). (180) 


But, since 


and 


dfl dll . . ri , 

r^=“C 0 si (Equation 51) 

fit V 

dE = n^a 0 A> (Equation 1?4) 


we can write Equation 180, after multiplying 
through by dt./dE, as 

dX, x i 1 h . | . r dSl .... _ s , 

dE^ 2+ 2n 0 (1+I,) r d^ C0S ^^-mX a )A. 


But we know 


%=a 0 3 ' 2 , 


and 


A 0 


1 

A a 0 (l-Co 2 )’ 


r=r(l + v), 


so finally we have 

d\i . , 1 h a 0 dQ ....... . .. 

-7y-,=aX 2 +o r ~= - — N 7 cos i(fx 4 — mX 3 )A. 
dE 2 A„ v i_ fo 2 dE 

(181) 


The derivatives of the other three X parameters 
are developed in precisely the same manner as 
was d\i/dE. The final forms they assume are: 


mined in the previous iteration for these X’s. 
Also, it is readily seen that every component of 
the derivatives can be expressed as a trigonometric 
series in E and (to). No F’s will occur anywhere 
in these derivatives. These are the final forms 
which are formally integrated to obtain the new 
values of the X’s in any iteration. 

EXPRESSIONS FOR E, AND v 

h h o 

We have shown previously (Equations 135 and 
137) that the IT function can be separated as 
follows: 

IT=E+T cos E’+'E sin F, 

where E, T, and T are series containing E and (to) 
alone, and are thus independent of F. And we 
have also shown that 


3 +- T — '-r+r 


If we set 


- = 1+A, 


(186) 


we can write 


f = 1-A+A 2 -A 3 + . . 
n 0 


(187) 


which gives 

E+e„T=-l-(l + A) + 2(l-A + A 2 -A 3 + . . .), 

or 


dx ; 

dE 


= — aXi + 


1 h a 0 dtt 


2 v l-e 0 2 &P 
i u dll 


dX 3 .1 h a 


cos i(— X 3 l— X 4 m)A, 
(182) 

cos i (X 2 f+Xim)A, 

(183) 


E+e 0 T=-3A+2(A 2 — A 3 + A 4 - . . .). 
Writing this in a form suitable for iteration we have 

A=-| (S+c„T)+| (A 2 — A 3 TA 4 — . . .). (188) 


r/X 4 . , 1 h a 0 Off v x ; i x \ . 

m = ~^ + 2 k ,, ,. o i d^ cos t(_XiZ+X2m)A - 

(184) 

In each of these four derivatives, cos i is expressed 
as 

cos f=X 4 2 +X 3 2 — X 2 2 — X 4 2 . (185) 

By inspection, we notice that these four expres- 
sions are set up in such a way that they allow 
iteration. Each derivative is expressed in terms 
of the X’s themselves, and uses the series deter- 


From Equations 136 and 186, it immediately 
follows that 

^=1+2 (*+ f oT+A). (189) 

The series for IT determines E and T. Then, 
from Equations 186, 188, and 189 the remaining 
quantities A, h 0 /h, and h/h 0 are found. This 
method of iteration is desirable in that it avoids 
the use of further integration to determine these 
series. Once the series h 0 fh and h/h 0 are found, 
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the perturbations of the radius vector v can be 
readily found from our original definition of W, 

IX (Equ “ ion 125) 

which, after substitution of A 0 /A= (1 + A) can be 
written suitably for iteration as follows: 

r=i(A-F)-^(lF+A>. (190) 

Again we are avoiding additional integration by 
the use of this iteration process. The W series 
used in this iteration process is obtained simply 


by replacing the F’s by E’s in the series generated 
for W. 

We now have developed all the necessary 
expressions which are used in the theory. In 
Equations 173, 154, 181, 182, 183, 184, 186, 189 
and 190, we have developed series forms for 

dW dn 0 Sz d\i d \ 2 d \ 3 rfX 4 h 0 h , 
dE’ dE’ dE’ dE’dE’h’ ho and 

respectively. Our next problem is to add the 
constants of integration which are determined in 
very specific ways; the method of determination 
of these constants has led to much of the contro- 
versy and confusion circulating about Hansen’s 
theory. 


SECTION VII 


DETERMINATION OF THE SECULAR MOTIONS y, a, AND r, 
AND THE CONSTANTS OF INTEGRATION 


DISCUSSION 

In Musen’s original paper (Reference 6) he 
points out that the “real constants of integration 
are the six elements a 0 ,e 0 , g 0 , do, and co 0 =x 0 — a 0 .” 
It is understood that any theory which attempts 
to give the motion of an orbiting body in a disturb- 
ing force field requires a complete solution of the 
differential equations for the time rate of change 
of the six elements which define the orbit. The 
constants of integration of the set of six differential 
equations would clearly be some starting values 
of the six elements. In Hansen’s theory, the 
differential equations of the elements are disguised 
and combined in the differential equations for 
W, n 0 dz, and the X’s. Intuitively, therefore, we 
would expect the constants of integration of these 
new differential equations to be functions of the 
elements alone. This is another way of stating 
that Hansen’s theory introduces no new constants 
of integration, other than the elements. 

Musen’s next statement is that these elements 
“do not have any simple kinematical or geo- 
metrical meaning; in particular, no moment of 
time exists for which these elements are osculat- 
ing.” The point here is that the “real constants 
of integration” are mean values of the osculating 
elements, which Hansen assumes to be known ex- 
actly. However, in practice, the mean values are 
never known as an initial condition. 

The standard technique for determining them 
is to assume a set of elements at some time t Q , use 


them to predict an orbit, and then to make an orbit 
correction assuming that residuals are a result of 
inaccurate starting values of the elements. Then 
by some correction device, new starting values of 
the elements are found and the process is repeated. 
In general, the first approximation is the set of 
osculating elements determined by observation at 
some time f 0 . After several corrections, the values 
of the elements should converge to the “mean ele- 
ments” which are Musen’s “real constants of in- 
tegration.” However, since these elements do not 
appear explicitly in the constants of integration of 
our development, let us turn to the matter at hand. 

At this stage of the exposition, it behooves us 
to recall the basic ideas that have governed the 
development thus far. After having introduced 
an auxiliary ellipse into the osculating orbit plane, 
we established that the ellipse moves uniformly in 
this XY plane with respect to the eccentric anom- 
aly. That is, the perigee moves with an angular 
velocity y(dE/dt ) in the XY plane. The pertur- 
bations of the real satellite’s motion in the orbit 
plane are combined in the two perturbations, v 
and n 0 Sz. The perturbations of the orbit plane, 
that is, of the elements d, <r, and i, are contained 
in the four interdependent X parameters. By the 
way in which we introduced and defined the secu- 
lar motions y, a, and rj, we have established the 
following important condition: The four param- 
eters Xi, X 2 , X.-j, X 4 and the two perturbations v and 
nodz can contain no secular terms. The secular 
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motions a and 77 are defined as containing together 
all the secular motions of <7 and d, and we know 
that i physically has no secular motion. Also, v 
clearly has no secular motion, and the secular 
motion of nobz is all contained in y. 

Since the development uses the iteration process, 
the above condition dictates the evaluation of the 
various constants of integration. In most cases, 
if one of the above six series contained a constant 
term, integration of the series would produce a 
secular term, which is not allowed. This is the 
logic behind the determination of the secular 
motions, and the arbitrary constants of integra- 
tion. 

DETERMINATION OF y 

Now let us examine each integral separately. 
The first integration performed is that of d\V/dE. 
As was explained in Sections IV and V, we finally 
were able to develop the differential equation 
dW/dE as a trigonometric series whose arguments 
are of the form [iE-\-j(w) f-kF\. From the inte- 
gration of this derivative, we will be able to per- 
form the important determination of y, the secular 
motion of the perigee. The condition that allows 
this result has been given: nobz cannot contain 
secular terms. 

The determination of the constant of integra- 
tion of dWIdE is a result of the condition stated 
above, as well as the additional condition that 
n 0 bs can contain no term of the form sin E. This 
second condition derives simply from t lie original 
definition that n u bs is the deviation of the real -satel- 
lite from Kepler’s equation, and Kepler’s equation 
already contains a sin E term. It should be 
emphasized that this is strictly a result of the way 
in which nobz is defined. We will presently show 
that these two conditions dictate that the constant 
of integration in IF be Co+Ci cos F. 

Thus, the final form of H" is as follows: 

cos [iE+2j(u)+kF] 

i j k 

s.k s i n [iE-]- (2J+ 1 ) (co) +&/'’| 

i J k 

+C0+C1 cos F, (191) 

where k= — 1, 0, +1. 

The limits on k and the form of the coefficients 
of (co) are not arbitrary but are purely results of 
the way in which various terms combine. The 
expression above is the general form of IF in 
every iteration. We need to know this form 
before we can explain the determination of y. 


Now, as can be seen above, and in accordance 
with the condition that n„bz have no secular terms, 
IF is not allowed to contain any secular terms. 
Therefore, since the integration is performed with 
respect to E, all terms which have F alone in the 
argument must be removed before integration. If 
they were not removed, W would contain terms 
of the form (sin F)E, a secular term. There are 
no constant terms in dW/dE, as will be seen by 
inspection of Equation 173, and so no secular 
terms arise from this source. We shall presently 
show that the only term of the series for dW/dE 
which contains F alone in the argument is of the 
form {Ax + yA-f sin F where A, and A> are con- 
stants. Thus, y is determined in such a way 
that this term disappears. 

To show that dW/dE contains a term of this 
form as the only term in F alone, we must examine 
Equation 173 which is 


dW 

dE 


■■NAr 


da n U 

nr 


F.1/A 


d a-pQ Sy 
dE fl — e 0 2 


We know the following: 

1. The quantities r(da f) Q/dr) and 5 a 0 QJdE are 
series containing arguments of the form 
iE-\-j(u>) only, because both are the results 
of “barring” r(b£l*/dr) and dtt*/dE, that is, 
replacing F by E. 

2. Also, IF and v are series in E and (u>) alone, 
since IF is IF with the F’s replaced by E’s 
and v is a function of IF. 

3. The quantity r la n =l—e 0 cos E. 

4. The quantity p/a 0 = 1— e 0 cos F. 

5. The quantities h u /h, h 2 jhf, and 1/(1 -F n) are 
series in E and («) alone since they arc all 
functions of IF. 


We will now show that the products NAr(da 0 tt/dr, 
and MA(da n ( .l!dE) can have only terms in ±F 
alone, that is, no terms of ±kF can appear with 
|/t|>l. After this we show that NAr(da 0 tt/dr) 
and A/A (daAlldE) can contain only sin F terms 
and can have no cos F terms. 

Examining NAr(da 0 tt/dr) , we can see from items 
2 and 3 above that A contains only terms of E 
and (co) since 

a=±=£A+^=^ 

1+H \ x 1 -c 0 2 11 0 
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and we know from item 1 that r ( da All dr) contains only 
E and (co), not F, Furthermore, we know that M 
and N contain terms of the form cos (mE-\-nF) 
where — and n= — 1 , 0 , 1 (see Equa- 

tion 173). This means no multiples of F other than 
— 1 , 0 , 1 exist in terms of M and N. Therefore, 
it is evident that the product A 7 A r(da 0 Q/dr) contains 
no terms with multiples of F other than — 1 , 0 , 1 . 
Exactly the same arguments can be used in the 
case of MA(daAl/dE) . 

Having seen that no terms of arguments kF 
with |*|>1 can appear in the first two terms of 
dW/dE, we want to show that only sin F terms can 
appear, and no cos F terms are possible. The 
proof of this rests squarely on two facts that are 
ascertained by close inspection of the series forms 
that go into ft and A. This inspection shows the 
forms to be as follows. (In these series expres- 
sions, (7 and S represent only general coefficients.) 

C i , J cos[iE+2j(u)] +X1S St. j sin 

i j i j 

[fE'-f (27 + I) («)] (See Equation 118). 

The two partials of ft then are easily found to be of 
the form: 


are produced: 

Ar Ar ^=^^ C 7 i , J .,cosff£+(2j+l)(a,)+^] 

Si.j.n sin [iFj-\-2j(w)-\-kF] 

i j k 

and 

A/A ww=2jZ1Z1 Ci.j.k cos [i,E-\-(2j-\- l)(co)+&E] 

Oii i j k 

Si.j.ts in [iE-\-2j(w)-\-kF]. 

i j k 

It is evident from these series forms that neither 
of the two terms can contain cos F terms since 
all cosine terms must also have an (w). The sine 
terms can clearly be in F alone, but only ±F, as 
we have shown previously. Thus, we know that 
only one term in the series for 

, ri da 0 ft , da 0 ft 
NAr +MA 

or oE 

can contain F alone, and that term will be of the 
form At sin F. 

Let us now examine the third term of dW/dE, 
that is, 


Ci.jC os [iE+(2j+l)(w)] 

-SS St,] sin [iE+2j(a)] 

i 

and 

r ™ Oi,, cos [iE+2j(o>)} 

+S2 Si, } sin [iE-\- ( 2 j — |- 1 )(<*>)], 

> j 

and the series A can be expressed as 

A=22j Ci,} cos [tE-\-2j(,w) 1 
< i 

+ 2 S St,] sin [iS+( 2 j+l)(ti>)]. 

i j 

We know further from Equation 173 that M and 
N take the general form 

M= y~)T" ) Ci,k cos ( iE-\-kF ) 


— where S= p W -|- 1 + 7 ^^) e 0 sin F. 

yl — e 0 2 «o 5/ V hj 

(See Equation 173). 

We can examine S alone, since y and e 0 are con- 
stants. First let us look at IF. We can show, 
by examination of Equation 191 that 

IF=Co + Ci cos F -\- Periodic Terms, 

but contains no terms other than c! cos V which 
have F alone in the argument. We know this 
because all such terms are removed from dW/dE 
before it is integrated . The Ci cos F term is part 
of the added constant of integration. From this 
fact we know that 


c i s i n E-)- Periodic Terms, 

again containing no terms in F alone other than 
—Ci sin F. So we see that, since 


and 

N= y^.y^, S t 1 sin (iEAkF). 

i k 

Using these five general forms, we can see that 
the following general forms of the product series 


— =1 — c 0 cos F, 


the first term of S, 


1<A¥ 

a 0 dF’ 
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can be written 

p dfF ,, . ,, 

— ir™=(l — «o cos F)(— C! sm f 
do Or 

+ Periodic Terms in E and F) 

= —c i sin Cie 0 sin 2 F 

-(-Periodic Terms in E and F. 
The second term of S is 

+ e 0 sin F=We 0 sin F-\-e 0 sin F 

+j e 0 sin F 

=c 0 e 0 sin F-\-Cie 0 cos F sin F 

-f e 0 si n E-j- Periodic Terms in E and F, 

since h 0 /h is a series in E. Therefore, S becomes 

$=[ — ci— e 0 {c 0 +l)] sin F 

+ Periodic Terms in E and F. 

So, 

dlV 

-jg=(A l +yA i ) sin F 

-(-Periodic Terms in E and F, (192) 

where 

~ ^2=[c 1 +e 0 (c 0 +l)]Vl— e 0 2 . 

Thus, from Equation 192 we are able to see that 
proper adjustment of y leaves dWJdE with the 
desired periodic terms in E, (co), and F. We now 
turn to the constant of integration problem. 

CONSTANT OF INTEGRATION OF W 

As we have stated before, the condition that 
dictates the form of the constant of integration in 
W is that n„i >2 is defined to be the deviation from 
Kepler’s equation: 

E— e 0 sin E=g 0 +n 0 (f— t 0 )+n 0 S 2 . (Equation 75) 

Therefore, since a constant term, a secular term, 
and a term of the form sin E all appear in Kepler’s 
equation, none of these is allowed in n 0 Sz. 

Returning to the differential equation (Equa- 
tion 154) we have 


dn 0 Sz 

IF+r 2 r I 

a-r 2 \ y r 

r \ 2 

dE 


vl+B7Vl-c 0 2 V 

ij 


which after minor rearranging can be written 

dn 0 8z _yy r | f v 2 —W 2 f_ y l—v 2 r 2 ^ 
dE ~ do y i +w a 0 yjAT^A 1+ w a 2 ) 

(193) 


From this equation, we can derive a form con- 
venient for the use of the iteration process. In 
any given iteration loop, the factors and the terms 
in parentheses are all of the previous iteration, 
with the exception of y which is the latest value 
found. For IF in the first term on the right-hand 
side of Equation 193 we "bar” operate on Equa- 
tion 191 and get: 

W=c 0 +Ci cos E+Y^^Cij cos [iE+2j(w)\ 

+EE«u s in [iE+(2j+l)(«)]. (194) 

» i 

It should now be clear that we set the form of 
the constant of integration in W to be c 0 +C! cos F 
in order that we could use it to remove the con- 
stant terms and the cos E term of the series for 
dn 0 8z/dE. If not removed, these terms would 
give secular terms and a sin E term after integra- 
tion. Thus, c 0 and c, are set in such a way as to 
make the terms disappear before the integration. 

The subsequent expressions for c 0 and c i are 
obtained as follows. Since we are using previously 
generated series for the terms in parentheses in 
Equation 193, and the new series for W in the 
term W(r/a 0 ), we have: 

^pw==(co+Ci cos ^)( 1— e o cos E): 

+S Z) G i. i cos HE+2j(w)] 

i j 

+2 2 S t . , sin \iE+ (2 j+ 1) («) ], (195) 

i j 

since r/a 0 =l— e 0 cos E. By means of the stand- 
ard trigonometric identities, Equation 195 becomes 

<n«o+(ci — C 0 c 0 ) cos E — i c,c 0 COS 2 E 
+2 S C { .j cos [iE+2j(<a)] 

i j 

+Z s i n [iK+(2j’+l)(co)]. (196) 

i j 

The two conditions that the constant terms and 
terms of the form cos E must be eliminated dictate 
that 

c 0 -|c 1 e 0 +C'o,o=0 (197) 


and 


Ci — CoCod - ( V o — 0- 


(198) 
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And a further consequence is that the coefficient 
of cos 2 E in Equation 196 must be added to the 
coefficient 0 2 , o of the series. From Equations 197 
and 198, we are able to write 


out again that in Equations 201 and 202, the X’s 
that appear on the right-hand sides are X’s of the 
previous iteration.) 

If we let 


c 0 =-(^|£t¥^ 0 ) (199) 

and 

C i = _(2co^±^Lo). (200) 

[Note: These values can also be found by itera. 
tion. See Appendix B, steps 61-67.] 

Putting these values into the series (Equation 
194) gives a final series form of IP. We now are 
able to write that dn 0 8z/dE consists of periodic 
terms only. When dn 0 5z/dE is integrated, the con- 
stant of integration in n 0 Sz is set equal to zero, as 
a result of the condition that n 0 8z contain no con- 
stants. And so we are able to develop series for 
W and n 0 8z for each iteration. We now must 
turn to the secular motions a. and 77 and the con- 
stants involved in the integration of the X differ- 
ential equations. 

DETERMINATION OF a AND v 

The differential equations for the X’s, Equations 
181 through 184, provide an easy determination of 
a and 77. Recalling that the X’s were defined in 
such a way that they contain no secular terms, it 
is at once clear that the derivatives can have no 
constant terms. The expressions for d\ 2 !dE and 
d\ 3 /dE are especially convenient here. They are 
given by 


dX 2__ .1 A a 0 bfl 

dE 1 2 h 0 ^J\ — e 0 2 


cos i(— \ 3 l— \ t m)A 


( 201 ) 


and 


tfX 3 . . 1 h a 0 dll 

dE~ + vXi+ 2h 0 ^TI^di 


cos i(X 2 £+Xim)A. 

( 202 ) 


We know that the series for X 4 and X 4 contain con- 
stant terms, since the first approximations to Xj 
and X 4 are sin (i 0 /2) and cos (i 0 / 2), respectively. 
We also know, from the way they are defined, that 
X 2 and X 3 can contain no constant terms. The 
first approximations are X 2 = X 3 =0. (We point 


rp 1 fb &o 50 / .s 

1 =0 T -7 XT ( cos V A , 

2 A)Vl-c 0 2 

(203) 

&i=constant part of Xi, 

and 

(204) 

/; 2 = constant part of T(— \ 3 l— X 4 m), 

, (205) 

then we remove the constant part of dX ; 

ifdE by 

setting 

r 


02 

a= T 

(206) 

In a similar manner, taking the expression for 

d\ 3 [dE, Equation 202, and letting 


& 3 =constant part of T^f+X!?/!) 

and 

(207) 

h 4 = constant part of X 4 , 

(208) 

we have 


V 64' 

(209) 


We now are able to integrate formally all the 
X derivatives. Upon integration, we are able to 
set the constants of integration in X 2 and X 3 equal 
to zero, since our conditions dictate that X 2 and 
X 3 contain periodic terms only. It is a little more 
involved to find the arbitrary constants in X 4 and 
X 4 . Two further conditions must be satisfied, and 
these govern the arbitrary constants in X 4 and X 4 . 
They are: 

1. That X! 2 +X 2 2 +X 3 2 +X 4 2 =l. (See Equation 
89.) 

2. That the principal term in the latitude must 
have the form sin i 0 sin [/+(«)]. (See 
Equation 87). 

Since we had the form of the latitude 
2(X 4 X 4 — X 2 X 3 ) sin [/+(&>)] 

+ (2X 2 X 4 + XA 3 ) cos [/+(“)], (see Equation 91) 

condition 2 dictates that the constant part of 
2(X 4 X 4 — X 2 X 3 ) be sm %$. 

Now, we have stated previously that parts of 
the constant terms of Xi and X 4 are the first approx- 
imationss, in (i 0 /2) and cos (i 0 /2), respectively, 
and we know that neither X 4 nor X 4 can contain 
secular terms. Let us arbitrarily pick a form for 
the constants of integration which will enable us 



38 


TECHNICAL REPORT R-14 7 NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


to satisfy the two conditions above, 
the final forms of Xi and X 4 : 

We write 

Xi^sin -y+q (A-\-B)~ | 5X 1; 

(210) 

X 4 =cos {A— 5)+6X 4 , 

(211) 


where A and B must be determined ; and 
cos [ 'iE+2j(w)] 

i j 

3 si n [i-E+(2J+l)(w)], (212) 

< j 

5X 4 =SSC , i.j cos \iE+2j(w)\ 

i 3 

+ 'T','F', St i sin [i£'+ (2j+l)(w)]. (213) 

i j 

These forms of Xi and X 4 are actual series forms 
arrived at after several iterations. Likewise, we 
find that 

*X 4 =£!S&.j8in [iE+2j(w)\ 

+2SQj cos [^£'+(2J+1)(co)], (214) 

i j 

SXs^SS'S'm s >>‘ [i-E’+2j(w)] 

» i 

i cos \iE-\- (2/+ 1 )(&>)], (215) 

i i 

where in final form, 

X 2 =5X 2 , (216) 

X 3 =5X 3 . (217) 

Taking the final forms of the X’s (Equations 210, 
211, 216, and 217), we apply them into the two 
conditions 

Xi 2 TX2 2 TX 3 2 -(-X 4 2 = 1 (218) 

and 

constant part of 2(X 4 X 4 — X 2 X 3 ) = sin i 0 . (219) 

In this procedure we must remember that even 
though X 2 and X 3 do not contain constant terms, 
their product terms will. The following results 
are obtained by setting the constant parts of 
each side of Equations 218 and 219 equal, and 
by setting the periodic parts of each side equal. 
Clearly, in this manner, the periodic parts vanish 
from the equations and we have the forms from 
which A/2 and Bj2 are determined. 


After making the substitutions, we find that 

2 (A iJ r IF) -\- A ^cos ^+sin B ^cos sin 

+ [constant part of (5 Xi 2 +5X 2 2 +5X 3 2 +5X 4 2 )] = 0 

( 220 ) 

and 

2 (A 2 — B 2 )+A (cos ^+sin ^)+B (cos sin 

-(-[constant part of 2(5Xi(iX 4 — 5X 2 5X 3 )] = 0. (221) 

Solving Equations 220 and 221 for A and B, we 
get: 

A 2 +2A ^cos -^T-sin —^-(-constant part of 

[(5X l +5X 4 ) 2 +(5X 2 — 5X 3 ) 2 ]=0 (222) 

and 

B 2 —2B (cos sin -^)+constant part of 

[(5X 1 -5X 4 ) 2 +(5X,+5X 3 ) 2 ] = 0. (223) 

It is preferable to solve for A and B by iteration 
rather than by t lie quadratic method, so we re- 
write Equations 222 and 223 as 



where H and G are the bracketed terms of Equa- 
tions 222 and 223, respectively, and the A/2 and 
B/2 used in the right-hand sides are the values ob- 
tained in the previous iterations. 

We have now described completely the tech- 
nique and reasoning used in the determination of 
the secular motions y, a, and 17, and the constants 
of integration in W, no 8z, and the four X’s. These 
six formal integrations are the only ones which are 
performed. All other determinations use the iter- 
ation process, as we have shown. The next sec- 
tion attempts to outline in detail the entire pro- 
cedure used in the generation of the final series 
forms which are used in the eventual predictions. 
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SECTION VIII 


PROCEDURE USED IN THE GENERATION OF FINAL SERIES FORMS 

DISCUSSION 


Iii the preceding sections we have developed 
explicit expressions for all the equations which 
must be solved in the problem. It is now possible 
to discuss the entire procedure for the generation 
of final series forms, using the various expressions 
found. 

It is important to emphasize one point. 
Throughout the process of iteration, our goal is to 
ascertain final, “good” series expressions for the 
following functions: n 0 8z, v, X 4 , X 2 , X 3 , and X 4 . 
We start out with approximations to these series 
as follows: 

Jlo52 = 0, 
y=0, 

\ • 2*o 

X!=Sln -=> 

*2 = 0 , 

X 3 =0, 

X 4 = cos A. 

Out of the iteration process, we eventually 
converge to series which express accurately these 
six functions as functions of the eccentric anomaly 
of the fictitious satellite E and the mean argument 
of perigee (to), which is in turn a function of E. 
Out of the process come the three secular motions 
y, a, and 17 . Conceptually speaking, the entire 
procedure is one of separating all secular motions 
from the periodic motions. Thus, series are 
generated which give the motion of the fictitious 
satellite as defined, and give the relationship 
between the two coordinate systems. And so, 
the final value of y and the final series form of 
71 adz provide that all the periodic angular perturba- 
tions of the satellite in the orbit plane are con- 
tained in the mean anomaly, and the secular 
motion in the plane due to disturbing forces is 
contained in the secular motion of the perigee. 

We must remember, therefore, a fact that is 
often overlooked and consequently is the cause of 
consternation in studying this theory. None, 
absolutely none, of the series is evaluated until the 
final series forms are determined after several 
iterations. Only the coefficients of the trigono- 
metric terms are calculated; the terms themselves 
are not, but are carried along throughout the itera- 
tion process. The iteration process is used to 


arrive at final series forms expressed with numeri- 
cal coefficients multiplying trigonometric terms of 
indexed arguments. Once the final series are 
developed, we go into the final stage where we can 
predict the position and velocity vectors of the 
real satellite at some future time t. Thus, there 
are two major phases of the development: the 
generation of final series forms, and the calculation 
of position and velocity at various times. 

The first of these phases has been written as a 
machine program by G. E. Collins of IBM, in col- 
laboration with Paul Herget of the Cincinnati 
Observatory. Its IBM code name is the General 
Oblateness Perturbations (GOP) program. The 
second phase, considerably more straightforward 
than the first, has been written as a program by the 
IBM Space Computing Center and is known as the 
IKIXT program. Now, let us trace through the 
procedure used in the first of these two phases. 

For any given satellite, we must begin with 
nominal values of the three elements a, e, and i. 
The theory calls for “mean” values of these oscu- 
lating elements, which we will assume we have. 
In actuality, these are determined by starting with 
nominal values, predicting an orbit, running an 
orbit correction, and then correcting the values of 
a, e, and i used. But let us say we have good nu- 
merical “mean” values of a, e, and i which we 
designate as a 0 , e 0 , and y,. Oidy these three ele- 
ments are needed in the development of the series. 
We also need values of the three geodetic param- 
eters ko, I- 3 , and & 4 . Next, we need the following 
first series approximations, each of which is used 
in an iteration process: 


Xi=sm 

O 

II 

<0 

0 

X 2 — X 3 =0; 

ho . 
A 

X 4 =cos 

A_i 

W= 0; 

ho 

A~0 


and we need first approximations to the constants, 
n 1 A B 

V= 0 an d -== 2 = 0 . 
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The first step is to develop the series for dW/dE. 
This involves the partials dil*/dF and r(dQ*/dr) 
which are readily found, once Q* is developed in 
terms of \ p*, where the A’s are those listed in 
Equation 226 for the first approximation, after 
which the last series obtained for the A’s are used. 
Taking Equation 173, 


dW \ T * da„fl , daoO , Sy 


we see that N and M contain h/ho and l/(l + r), so 
we must use the series forms obtained in the previ- 
ous iteration. The same is true of the v, IT, and 
h 0 /h that appear in A and S. And finally, we see 
that y appears twice in dWjdE; as a coefficient of 
S' in Equation 173, and in the expression for A. In 
the expression for A, the value of y previously 
determined is used, whereas the y multiplying S 
is the y to be determined by the removal of sin F 
terms from dW/dE (see Section VII). 

Once dWjdE is determined and the sin F terms 
are removed, it is formally integrated and the con- 
stant of integration coT<T cos F is added. Keep- 
ing Co and Ci undetennined_as yet, we “bar” the 
new series for IT to get IT which has now the 
terms c 0 and c, cos E. Keeping this series for W 
available, with c 0 and Ci still undetermined, we 
turn to our expression for dn^z/dE 


dnd&z jp r / v 2 — IT 2 7 y 1-r 2 P\ 
dE H «o + \l+W «o /l — e 0 * 1 I IT «oV 

(Equation 193) 

Into this expression we put the v series deter- 
mined in the previous iteration (or, in the case of 
the first iteration, e = 0) . In the term in paren- 

theses, we also put for IT the series determined 
previously. The y is the number just found prior 
to the integration of dW/dE. Then, in the first 
term, IT(r/a 0 ), we put the IT series just deter- 
mined and find co and Cj so that no constant terms 
or terms of the form cos E appear in dn 0 8z/dE. 
The expressions for c 0 and Ci are given in Equa- 
tions 199 and 200. With Co and Ci evaluated, we 
formally integrate to get the series for n 0 Ss, which 
contains no additive constant of integration for 
reasons given in Section VII. We now also have 
the new series for IT and W since c 0 and Ci have 
been determined. 


Our next step is to take the new series for IT and 
write it in the form, 

IV=E+T cos E+'T sin F. 

This is accomplished by scanning the series IT 
and first picking out all the terms which have no F 
in the argument, a series we call E. Next we are 
able to find T by putting F= 0 in all the terms in 
which F occurs. It should be remembered that k 
only takes on values of — 1, 0, 1 in the series for IT 
(Equation 191). This splitting process simply uses 
the trigonometric identity for the sine or cosine of 
the sum of two angles. With the series for S and 
T determined, we turn to our equation for A 

A =-| (H+c 0 T)+| (A?-A 3 + . . .). 

(Equation 188) 

Using the value of A of the previous iteration, or 
for the first iteration A=0, in the right-hand side, 
we solve to get a new series for A. Properly, for 
any iteration, the expression should be written 

A „+t=-5 (E+e„T)+? (A„ 2 -A„ 3 + . . .). (227) 

In any given iteration of the whole program, that 
is, when new series for IT, IT, n o5,r, etc., arc found, 
we iterate around Equation 227 n times, until 
|A b+1 — A„[<(e where e is some arbitrarily chosen 
small number. This A„ +1 is the A to be used in the 
determination of h/ho, h 0 /h, and v. 

With the A series determined, and with the E 
and T series used above, ho/h and h/ho are found 
simply by substitution in 

^=1+A (Equation 186) 

and 

yr -= l-| _ Q (A-j-S-t-CoT). (Equation 189) 

These series for ho/h, h/ho, and A are then put aside 
with IT and n<>8z for use in the next overall iteration. 
The last series to be found is that for v, 

(A— IT) - x (ITT A) in (Equation 190) 

We use the same type of iteration of this equation 
that was used to find A. Here, the A used is the 
one just found and stored, and IT is simply the 
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last W determined. Again, the expression should 
properly be written 

r„ +1 =! (A-W)-| (W+A)p„ (228) 


grated. No constants of integration are added 
to X 2 or X 3 , for reasons previously discussed, so 
the integrated forms are stored as the new X 2 
and X 3 . 

The constants of integration added to 


where n is the number of the iteration of the 
above equation. The series v n originally put on 
the right-hand side is that determined in the last 
iteration of the whole program. 

THE X DERIVATIVES 

We can now turn to the computation of the 
X series and their constants of integration. Taking 
the expression for dXi/dE as an example, we can 
illustrate the procedure used in all the four 
derivatives (Equations 181 through 184). Thus, 


d\ i_ 1 h Oo dO 


cos \ 3 m)A. (229) 


It is readily seen that the factor 


1 h 


a 0 


2 K .g Q 2 £h[/ 


TT (COS l)A 


(230) 


occurs in each of the four derivatives. This factor 
is developed using the h/h 0 series determined in 
the previous iteration, not the one found in the 
steps immediately before this. We use the previ- 
ous series to assure consistency through the overall 
iteration process. The dfi/di/' was found by oper- 
ating with the “bar” operator on dQ*/d\p* and 
cos i we have noted to be 


cos i — X 4 2 — |— X 3 2 — X- 2 " — X[ 2 , 

where the X’s are those resulting from the previous 
iteration. The A factor is the same as that used 
in the formation of dW/dE. As was true of the 
X’s that occurred in the expression for cos i, all 
the other X’s on the right-hand side of Equation 
229 are those of the previous iteration. This 
leaves only a undetermined. 

The determination of a and 77 was described in 
Section VII, on the basis that dX 2 /dE and d\ z jdE 
must contain no constant terms. In any itera- 
tion, a and 77 are determined as described before 
any of the X derivatives are integrated, since the 
expressions contain the value of a or 77 just deter- 
mined. Thus, dX 2 /dE and d\ 3 /dE are formed 
before d\,/dE or d\ 4 /dE. Once a and 77 are 
found, the four derivatives are formally inte- 



have also been discussed previously. The forms 
are 

Xi=sin 7^+2 (^1+-®)+ (231) 

X 4 =cos 7^+2 (-‘4"“-®)+ j* (232) 

where A and B are determined by iterations of 
Equations 224 and 225. As was described in 
Section VII, H and G in the latter are functions of 



C i /X, 
JdE 


before the constants of integration are added. 
The iteration process is to be handled in exactly 
the same manner as in the cases of A and v, so 
we could write Equations 224 and 225 as 



After the iteration process gives convergence to 
some A/2 and Bj 2, the new X x and X 4 series are 
given by Equations 231 and 232. 

We have now traced through one entire iteration 
of the first phase of the solution. We have gene- 
rated in this iteration new series for W, W, n 0 8z, 
hjh, h/h 0 , v, A, Xi, X 2 , X 3 , and X 4 . We have also 
evaluated new values of the constants y, a, r/, 
Aj 2 , and Bj 2 . With the iteration complete, the 
test for convergence is made with the secular 
motions. If 

\y n+i y n\ At, 


and 


«»+ 1— 

I’Jn+l — Vn\</e> 
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(where e is some arbitrary small number) the 
series are considered to have converged. If not. 
the entire process is repeated as outlined. This 
completes the description of the first phase, that 
incorporated in the IBM program known as 
GOP. Appendix B has been included as an aid 
to those interested in the programming of this 
development. It lists one possible method of 
attack but does not represent the only valid 


computational procedure. As can be seen from 
Appendix B, this program develops the X’s before 
it finds the series for h 0 /h,h/h 0 , and v. This order 
is preferable in that the h/h 0 series used in the X’s 
is that of the previous iteration, as discussed 
earlier. 

We next turn to the second phase of the solu- 
tion, that of producing a position vector of the 
real satellite at certain values of t. 


SECTION IX 


DETERMINATION OF FINAL POSITION AND VELOCITY VECTORS OF THE REAL 
SATELLITE AND ITS OSCULATING ELEMENTS 


INTRODUCTION 

The entire development thus far has led us to 
series which express exactly the perturbations of 
the real satellite and the motions of the orbit plane 
as functions of the eccentric anomaly. We now 
turn to the second phase, the determination of the 
position vector as a function of time. We know 
the motion of the fictitious satellite to be given by 

r sin f=a 0 -Jl — e 0 2 sin E 

and 

r cos /=a- 0 ( cos E—e 0 ). 

Also, we know the relationship between r and r to be 
r= (1 + r)r, 

and the relationship between the eccentric anomaly 
E of the fictitious satellite and the real time t of 
the real satellite to be in Kepler’s equation 

E—e 0 sin E=() 0 -\-n 0 (t— t 0 )-\-n<,8z. 

Thus, we are in a position to find the radius vector 
of the real satellite in the rotating XYZ coordinate 
system. However, we have to perform a rotation 
operation in order to express the position vector 
in terms of xyz coordinates, that is, inertial co- 
ordinates. 

THE ROTATION MATRIX 

The XYZ system has been rotated through 
three angles with respect to the xyz system. 
If we rotate the position vector through these 
same three angles with respect to the xyz system, 
the new components in the XYZ system are the 
same as those of the nonrotated position vector 
in the xyz system. Therefore, if we can find a 
rotation matrix operator which expresses the 
rotation of XYZ with respect to xyz, we can oper- 
ate on r with it to get r in the inertial system. 


One important thing to note here is that, for 
purposes of developing this matrix, we want to 
rotate our X axis from the x axis to the line in the 
plane of the orbit, from the origin to the perigee. 
This allows us to express r in the XYZ system in 
two components: r sin / and ? cos f, where/ 
is the true anomaly of the ellipse, that is, the 
angle between the position vector r and the line 
from origin to perigee. 

Therefore, to rotate from the xyz inertial, or 
equatorial system, into the XYZ system where 
the XY plane is the plane of the orbit, and the 
X axis goes through the perigee, we rot ate through 
three angles. Starting with the X, Y, and Z 
axes lying along the x, y, and z axes, we first 
rotate the XYZ system about Z, through the 
angle d. This leaves the X axis lying along the 
line of the nodes. Next we rotate the XY plane 
about X, through an angle i, putting the XY 
plane coincident with the orbit plane. Finally, 
we rotate the x and y axes about Z again, until 
X passes through the perigee. Clearly, this 
final rotation is through an angle (x„ +2/AK— a). 

A rotation about the Z axis through an angle a 
is given by 


A 3 [«]= 


'cos a —sin a 0" 
sin a cos a 0 

0 0 1 


(233) 


and a rotation about the X axis through an angle 
j8 is 


n o 


M$\= 


o 


cos 13 


0 ‘ 

—sin 13 . 


(234) 


0 sin /3 cos /3 
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Therefore, our total, or triple rotation, desig- 
nated by T, is clearly: 

T=A 3 [e]A[i\A 3 [ir 0 +yAE-a]- 

But from Equations (82), (85), (71), (83) and (86) 
we know that 

d— (0) — N-\-K and (ir 0 -\-yAE— a)= (u)-\-N-{-K, 


to express A 3 [K— 7V]^4i[i]^4 3 [i£-|-Ar] as one matrix 
in terms of Xi, X 2 , X 3 , and X 4 . We shall call this 
matrix 

Making the transformations 


f K-N=2a 

|^+iV=28 


J K=a + 8 

N= 8 — a 


so we write 


we have 




knowing that ^4 3 [ai + a 2 ]=.4 3 (a! 1 )Aj(a 2 ). We wish This is given explicitly by 




cos 2 a cos 2/3— sin 2a sin 2/3 cos i 
sin 2a cos 2/3 -f cos 2a sin 2/3 cos i 
sin 2/3 sin i 


— cos 2 a sin 2/8 — sin 2 a cos 28 cos i sin 2a sin i 
— sin 2a sin 28+cos 2a cos 2 8 cos i — cos2asini 
cos 2 8 sin i cos i 


which we write as 


X n = cos 2 a cos 28— sin 2a sin 28 cos i 



Xu 

X [2 

Xi3 

=^sin 2 (£)+cos 2 Q 

;)] 


X 2 i 

x 22 

X 2 3 

— r cos 2 ( 

i) 


_X 3 i 

k 32 

^33_ 

L V 



cos 2a cos 28 

sin 2a sin 28 


or 


Throughout the following development of the 
X m , „ terms, we use the definitions 


% % • 
X!=sin - cos N, X 3 =cos ^ sin K, 


and 


% 

X 2 =sin — sin N, X 4 =cos - cos K, 
Z z 


K=a+/3, 
iV=8 — a, 


and we have 


Xi=sin - cos (8— «), 


% 

X 2 =sin - sin (8— a), 


X 3 =cos - sin (8+a), 


X 4 =cos ^ cos (8+a). 


=cos‘ 


= 008 “ 


Xn=cos 2 (cos 2a cos 28— sin 2a sin 28) 

+sin 2 ( -0 (cos 2 a cos 28+ sin 2a sin 28) 

1 (Jj cos 2(a+8)+sin 2 cos 2(8— a) 

Q) [cos 2 (a+8)— sin 2 (a+8)] 

+sin 2 [cos 2 (8— a)— sin 2 (8— a)]. 

Then, employing X 1; . . ., X 4 , we have 
Xn = Xi 2 — X 2 2 — X 3 2 +X 4 2 . 

The same type of substitution and manipulation 
gives 

Xi 2 = — 2(X 3 X 4 +XiX 2 ), 

Xi 3 =2(X 1 X 3 — X 2 X 4 ), 

X 2 i = 2 (X 3 X 4 X,X 2 ) , 

X 22 — X 4 2 — X 3 2 — Xi' + X 2 2 , 

X 23 = — 2(XiX 4 +X 3 X 2 ), 


Now, we put 
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X31 = 2 (XA3 -f- X 2 X 4 ) , 

X32 — 2 (X4X1 X0X3) , 

X 33 = X 4 2 — )— X3 2 — Xf — X 2 2 . 
We now have our T matrix 



“x„ 

X12 

X13 

r=A 3 [(0)] 

x 2 i 

X22 

X23 


_X 3 i 

X32 

^ 33 . 


(235) 


where g 0 , the mean anomaly at the epoch, must 
be given. In the first iteration, we solve Equa- 
tion 237 by letting n 0 8z= 0. This gives us a 
value E 0 , which we immediately put into the 
series generated for n 0 8z: 

7i 0 5s=Fourier series in 

^JE+j^u 0 +(y+a—r])E^—j(y+a—r))E, 0 J> 

(238) 

which becomes 


which represents the transformation of coordi- 
nates from the osculating system to the equatorial 
system. Consequently, we will get the position 
vector in the equatorial system by the operation 


'a 0 (cos E—e 0 ) 


r=(l + r)r 


a 0 sin E-yJ 1 — e 0 2 ’ 


(236) 


n 0 8z = Fourier series in [iE’o+jcop]. 

Taking the n 0 8z, a numerical value, and putting 
it back into Equation 237, we solve to get the 
second approximation to E 0 , which we put back 
into Equation 238 for a second approximation to 
n a 8z. We repeat this iteration process until 


0 


| (£„) n+1 - (£„)„!<*, 


since a 0 (cos E—e n ) and a 0 sin E\ 1— e 0 2 are the 
components of r in the A’ and Y directions, 
respectively, considering A" to be. along the line 
from the origin (focus) to the perigee. Once 
the r, y, and z components of r are known, any 
other coordinates of the satellite can be readily 
found by elementary trigonometry. 

DETERMINATION OF r FOR A GIVEN t 

It is immediately evident from Equations 235 
and 236 that we can find r for any E if we know 
«o, 0o, and E 0 since we have series for v and the 
X’s in terms of E and (co), where 

(co) = Wo+(y+0! — 77 )/?— (y + a — v)E b . 

We need 6 0 in the T matrix, where 

kh[(0)] = ^l-3[0o — (a + i))E'+ (a + rfJEV]. 

We have already found the final values for y, a, 
and r). However, to get the components of the 
position vector at some time t, we must deter- 
mine the value E which corresponds to that time 
t. As well, to evaluate any of the series, we need 
to know the value E a which corresponds to t 0 . 

We proceed as follows: Taking Kepler’s equa- 
tion, we find the E 0 at t 0 by iteration. In this 
case, Kepler’s equation gives 

E o —e 0 sin E B =y 0 +n 0 (to— t 0 )+n 0 8z, (237) 


at which point we take E 0 to be known for t=t 0 . 
Next, we turn to some time t for which we wish 
to know the position of the real satellite. Turn- 
ing first to Kepler’s equation, in the first approx- 
imation we let n n 8z=0 and solve for E. Then, 
we put E and E a into the scries for n 0 8z and solve 
for a numerical value, which we then put into 
Kepler’s equation to solve for a second value of 
E. We continue in this iteration process until 

\E n +i — E 

Once we have determined the value of E 
corresponding to t, we are able to compute 
numerical values for Xj , X 2 , X 3 , X 4 , (1 + r), and can 
find (u>) and (0) as well as a 0 (cos E — e 0 ) and 
«oVT — e 0 2 sin E. Thus, we can easily put all 
of these values into Equation 236, and the com- 
ponents of the radius vector of the real satellite 
in the inertial coordinate system will emerge. 

DETERMINATION OF THE VELOCITY VECTOR 
AND OSCULATING ELEMENTS AT TIME t 

Bailie and Bryant have published (Reference 9) 
the method by which the velocity vector of the 
real satellite and the osculating elements at some 
time t are determined. The relationships follow 
directly and simply from the development dis- 
cussed here. 
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Recalling that the W function was split up to 
give 

TF=3+ T cos F+iIi sin F, (Equation 135) 
Bailie and Bryant define the two quantities 

^ e ° + x(^r 9 T ’ (239) 

and 

7 (240) 


Further, they use the components of the r 
rotation matrix written as 



-(Px) 

(Qx) 

(PxT 

r= 

(Py) 

(Qy) 

(Ry) 


.(Px) 

(Qx) 

(Px). 


(241) 


where (P x ), (Q x ), etc., are the components in the 
inertial system of unit vectors P, Q, and R. The 
unit vector P is in the osculating plane and is 
directed toward the perigee, R is normal to the 
plane, and Q=RXP. 

Bailie and Bryant give the components as: 


(Px)= + (x 4 2 — X 3 2 ) cos [(«) + ( 0 )] — 2 X 3 X 4 sin f (co) + (0)] 
+ (X 4 2 -X 2 2 ) cos [(co) — (0)] 

— 2 X ^2 sin [(«) — ( 0 )], 


(P v ) = + (X 4 2 — X 3 2 ) sin [ (co) + (0)] + 2X 3 X 4 cos [ (co) + (0)] 

— (Xi 2 — X 2 2 ) sin [(co) — (0)] 

— 2 X ^2 cos [(co)— (0)], 

(Pz)= +2(X 4 X 4 — X 2 X 3 ) sin (co) + 2(X 2 X 4 +XiX 3 ) cos (co), 

{Qx) = — (X 4 2 — X 3 2 ) sin [ (co) + (0) ] — 2X 3 X 4 cos [ (co) + (0)] 

— (Xi 2 — X 2 2 ) sin [(co) — (0)] 

— 2XJX;. cos [(co) — (0)], 

(Q u ) = + (X 4 2 — X/)cos [ (co) + (0) ] — 2X 3 X 4 sin [ (co) + (0)] 

— (Xi 2 — X 2 2 ) cos [(co)— (0)] 

+ 2XjX 2 sin [(co) — (0)], 

(Qz) — -f-2(XiX 4 — X 2 X 3 ) cos (co) — 2(X 2 X 4 -|-XxX 3 ) sin (co), 


(B x ) — _ h2(XiX 4 +X 2 X 3 ) sin (0) — 2(X 2 X 4 — XiX 3 ) cos (0), 
(#„)== — 2(X 1 X 4 +X 2 X 3 ) cos (0) — 2(X 2 X 4 — XjX 3 ) sin (0), 
(fi z ) = X 4 2 +X 3 2 — X 2 2 — Xj 2 . 

After some straightforward analysis the fol- 
lowing results are obtained giving the velocity 
vector v of the real satellite and the osculating 
elements in terms of /3, y, the components of P, Q, 
and R, h 0 jh, h/h 0 , and the X’s: 



sin E 

1 — e a cos Eh 0 2 

Vl — c 0 2 cos E h Vl— e 0 2 y ’ 
1 — e 0 cos E h 0 ~'~ 2 

0 


(242) 


a 


ffl 0 (l — e 0 2 ) A.Y 
1 — /3 2 — y 2 \ h/’ 


e = V(3 2 + 7 2 , 


tan i- 


tan co 


tan 0= 


V( F Z ) 2 +(Q Z ) 2 

X 4 2 +X 3 2 -X 2 2 -X! 2 ’ 

p(Pz)+vm 

m*)-y(Pz)’ 

(pjm-miP') 

\Px)m-{QxKP.) 


(243) 

(244) 

(245) 

(246) 

(247) 


And from the osculating value of the eccentric 
anomaly 


IT+e 

V 1— e 


e , E osl 

tan — 77 - 

e 2 


/T± 

V 1 - 


+ c 0 + E 
- tan — - 
e„ 2 


tan f 

7 V 1 — Co 2 


(248) 


the osculating mean anomaly M is found in 
Kepler’s equation: 

M=E osc —e sin E 0iC . (249) 

The components of the three unit vectors are 
easily arrived at, but are given explicitly by Bailie 
and Bryant, in terms of the four X parameters and 
the mean values (co) and (0). 
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SECTION X 

EVALUATION AND CONCLUSION 


This theory contains certain areas in which 
trouble may arise for particular values of the ele- 
ments. The three areas of difficulty in this theory 
are (1) small eccentricity, (2) large eccentricity, 
and (3) an angle of inclination in the neighborhood 
of the critical angle. 

SMALL ECCENTRICITY 

The first of these difficulties is not a weakness 
of the theory itself, but rather of the way in which 
it has been adapted for machine use. We recall 
that in the determination of the secular motion 
y, the technique involved was that of finding y 
such that no sin F terms appeared in the series 
for dW/dE. We solved an equation for y (see 
Equation 192): 

A.+yA^O, (250) 

where A, and A 2 are numerical values of coeffi- 
cients. Upon close inspection of the explicit 
expression for dWjdE, it becomes evident that 
this equation can also be written as 

eA[ -\-y(eA' 2 ) =0, 
or 


since the eccentricity appears explicitly in the 
coefficients which are summed to give A,: 

A l =(e n a l - 1 re 2 a 2 J re<£‘a i J reJa i J r ■ . .) sin F. (252) 

From the last term of dW/dE (see Equation 173), 
we see that 

Ai — eoA'z sin F. (253) 

The difficulty arises, not specifically from the 
appearance of e 0 in the denominator of Equation 
251 but from the fact that in Equation 252 the 
machine is adding a series of very small numbers 
during which process a large error can accumulate. 
When this sum is divided by a small number, 
e 0 A' 2 , the accumulated error is greatly magnified, 
and can easily exceed error limitations. An edu- 
cated guess as to the lower limit of the eccentricity 
would be on the order of 0.001. 

LARGE ECCENTRICITY 

The problem with large eccentricities is that the 
series for a 0 /p and a 0 /r converge very slowly, 


owing to the presence of the factor 2/V'l— e 0 2 : 

,-- 2 - = (h+fi cos A+/3 2 cos 2 FA- . . .\ 

P Vl-e 0 2 V2 / 

(Equation 108) 

It would be desirable to have different series ex- 
pansions for a 0 fr and a 0 fp for large eccentricities. 
It is impossible to find one series expansion which 
will give reasonably fast convergence for the entire 
range of eccentricity, O^eo^l- Also, the factor 
1/(1 — Co 2 ) appears in the development, for example, 
in the M and N factors of the derivative of W. 
The factor places a limitation on e 0 . The numeri- 
cal procedure described in this development is not 
satisfactory for eccentricities larger than approxi- 
mately 0.90. 

THE CRITICAL ANGLE OF INCLINATION 

It is a physical fact involving the particular 
oblateness of the earth, that at some angle of 
inclination in the region close to 63.4 degrees the 
forces disappear which gives rise to the secular 
motion of the perigee of a satellite. This is a 
result of oblateness symmetry with respect to a 
plane passed through the earth at this angle to the 
equator. The angle will vary slightly from 63.4 
degrees depending upon the geodetic parameters 
used in the potential function. But, nonetheless, 
there is some angle at which the secular motion 
of the perigee vanishes, or, in our notation, 

y+a— rj = 0. 

This leads to a problem primarily in the integration 
of the series dn u 8z/dE; the same problem is present 
in all the integrations, but in those of dWjdE, and 
d\/dE, the coefficients of the terms involved are 
very small and the difficulty is not readily ap- 
parent. We recall the series dn 0 8zjdE is of the 
general form 

i cos [%E+2j( a )] 

-j-Sij sin [iE+(2j+l)(co)]} . 

(Equation 195) 

If we take the term of the series where i=0 and 
7 = 0, we have the sin (co) term, where we know 

(o) = O) 0 + (y + a — rj) AE, 
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with A E=E—E 0 . So we must integrate with 
respect to E, that is, 


Jsin («)JE= J* sin [co 0 + (y + a — ri)E 

— (y+a—v)Eo]dE, 

which gives 

J - / \JT? cos (“) 

y+a—n 


It is evident that this term becomes meaningless 
when {y-\-a— rj)—>0, as it does at the critical angle. 
This renders the theory unsatisfactory in a very 
small region about the critical angle. 

We see that the term sin (u) is actually a very 
long period term in the region, and a constant at 
the critical angle. It might at first seem possible 
to separate the constant part of the term from 
its periodic part, but the periodic part would have 
to be considered secular and disallowed first order 
secular motions would result in the X parameters. 
This weakness in the theory seems to be inherent 
and unavoidable. 

In the other integrations, these long period terms 
are of a very small magnitude, and the limitations 
due to series truncation have thus far made it 
difficult to evaluate their significance. 

ACCURACY 

The theory, as stated in the introduction, is an 
exact one, which includes all orders of perturba- 
tions. The accuracy is that of the geodetic param- 
eters k 2 , k 3 , and & 4 in the potential iunction. 
However, in practice the accuracy is greatly 
affected by the truncation of series and the number 
of significant figures carried in the machine. The 
great number of series multiplications taxes the 
storage capacity of even the largest machines, so 
the question of series truncation is a serious one. 


[UNITS 

It has been found that the Vanguard system of 
units is quite satisfactory in this development, 
although by no means is it the only valid system. 
The theory can use any system so long as it is used 
consistently throughout. The Vanguard units 
allow that the product of the mass of the earth and 
the universal gravitational constant is unity. The 
unit of time, also a Vanguard unit, is the amount of 
time it takes an orbiting satellite at a distance of 
one mean equatorial radius from the earth’s center 


to travel one radian. This is 

1 Unit of Time=806.814 mean solar seconds. 

The unit of length, considered one mean equatorial 
radius, is 

1 Unit Length=6,378.165 kilometers. 
CONCLUSION 

The theory described has been adapted to make 
optimum use of the capacity and speed of modern 
computing machines. It is a purely numerical 
theory, the accuracy of which is determined only 
by the accuracy of the geodetic parameters, and in 
practicality, by the limitations of series truncation 
in the machine. The theory described deals only 
with general oblateness perturbations, but is 
easily adaptable to the inclusion of other per- 
turbing forces, including solar radiation pressure; 
solar and lunar perturbations, and the effects of the 
tesseral harmonics in the earth’s gravitational 
potential. Though this exposition has dealt only 
with the first three harmonics of the earth’s 
potential function, higher harmonics can be 
easily included. 
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APPENDIX A 


COROLLARY DERIVATIONS 


BASIC EQUATIONS OF AN ELLIPSE 

A basic property of an ellipse is 

a 2 =6 2 +c 2 , (Al) 

where a is the semimajor axis, b is the semiminor 
axis, and c is the distance from the center of the 
ellipse to one focus. The eccentricity of the ellipse 
is defined as the ratio 


e 


e 

a 


(A2) 


If we let the center of the ellipse be the center of 
a Cartesian coordinate system, and let the x axis 
lie along the semimajor axis, we can write the 
equation of the ellipse as 


x 

a 


2 

2 


y 

b 2 


(A3) 


and we can then solve for y, which is the distance 
r sin / (see Figure Al), by the equation 


(A4) 



Figure Al.— The geometry of the ellipse. 
We know as well that 

x=a cos E 

or 

x 2 =a 2 cos 2 E; 


y=^b 2 (l- J ~yrsmf. 
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therefore, we can write 

rsin j=Jb 2 ( 1- 


a 2 cos 2 E' 


=6V 1— cos 2 E 
= b sin E. 

But from Equations Al and A2, we have 


which gives 


b=^a 2 ( 1 — e 2 ), 


r sin/=aVl— e 2 sin E. (A6) 

Similarly, from Figure Al, we see that 

r cos/=a cos E—c; (A7) 

but since e=c/a, we can write this 

r cos/=a(cos E—e). (A8) 

Now, squaring Equations A6 and AS and adding 
them together, we have 

r 2 =a 2 [(l — e 2 ) sin 2 E+ cos 2 E— 2e cosi?+e 2 ] 

=a 2 [l+e 2 (l— sin 2 E)—2e cos E] 

— a 2 [ l—2e cos E-\-e 2 cos 2 E], 
so then 

r 2 =a 2 (l—e cos E) 2 

and 

r=a( 1—e cos E). (A9) 

Now, we can rewrite Equation A9 as 
1=- (1—e cos E). 

When this is multiplied through by e and re- 
arranged, we have 

fn 

e — — ( 1 — e cos E)= 0, (A10) 

which can be written 

e — cos E+cos E— ^(1 — e cos E)=0. (A10) 

But we know from Equation A9 that 


r 1 — e cos E’ 


so we can write Equation All as 

tp i ^ tp ^ ^ ? tp r (i , e a ' tp / * 

e— cos E -\ — cos E cos 2 E — cos E= 0. 

r r r r 

From Equation A8 we know that 
cos j=~ (cos E — e), 

and so we have 

e— cos cos E—^'j 

— e cos E (- cos E — — ^=0/ 
\r r / 

or 

e— cos E- fcos/— e cos E cos/=0. (A12) 

If we now subtract e from both sides of Equation 
A12, multiply through by e, and add 1 to both 
sides, we get 

1—e cos E+e cos/— e 2 cos i?cos/=l— e 2 , 

or 

(1—e cos E)( 1+e cos/) = l— e 2 . 
Multiplying through by a/( 1+e cos/), we have 

a(l— ecosE)=-^f i — (A13) 
1+e cos / 2 

but since the left-hand side of Equation A13 is 
r by Equation A9, we have the ellipse equation 

(A14) 

The next equation of interest states that the 
total area swept out per unit time by the radius 
vector of an ellipse is 1/h where h=ll\a(l —(F)- 

rdfi 


In polar coordinates, the area swept out in time 
dt is 

\ r(rdd ) 
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since the area is that of a triangle. So r 2 (dd/dt) 
is twice the area swept out per unit time. Now, 
the area swept out per unit time is the total area 
divided by the period of one revolution of the 
radius vector. Thus, 



(A 15) 


But A—wab where a and b are the semimajor and 
semiminor axes, respectively. We know from 
Equation A5 that b=Qa 2 { 1—e 2 ), so the area 
becomes 


are proportional to the sines of the corresponding 
opposite angles.” The sine of the geocentric 
latitude of a satellite, denoted by \ p, is defined as 
the sine of the angle between the earth’s equatorial 
plane and the line drawn from the center of the 
earth to the satellite. 

So, in the spherical triangle NQ8 (Figure A2), 
we have two sides (v— a) and SQ, and two angles 
i and SQN. If we designate the side SQ as a, the 
angle which subtends it, we can write the law of 
sines: 

sin SQN sin i 
sin {v—a) sin a 


A= TrWU-e 2 ). 

We know also that the period P is 2 ir/ra, where n 
is the mean motion (average angular velocity) 
of the radius vector, and from Kepler’s law, 
n=a ~ 3/2 . Therefore, P=2ir/a -3/2 and Equation 
A15 becomes 

A a ~ 3l2 =Qa{l-e 2 ) (A16) 

dt 2 7T 

which we have defined as 1/A. Thus, 

h= 1 = . 

Qa{l—e 2 ) 

Proof That ^=sin i sin (v—a) 

In spherical trigonometry, the law of sines is: 
“In any spherical triangle, the sines of the sides 


But we know that angle SQN is a right angle, 
and we have defined sin a as \p, so we have 

\p=sin i sin {v—a). 

DERIVATION OF EQUATIONS OF MOTION IN 
POLAR FORM 

Consider the motion of a particle in motion 
along a smooth path C, as shown. 




Osculating plane 
of the satellite. 


Distance VS = v 

Distance XN = <7 The arc SQ cuts the equatorial plane 

Distance NS = (v—a) perpendicular to it, by definition. 


Figure A2. — Hansen’s coordinate system. 




APPENDIX B 


THE COMPUTATIONAL PROCEDURE USED IN THE IBM GOP PROGRAM FOR THE GENERATION OF 

FINAL SERIES FORMS 


The following program was used in the actual 
computation of orbits. The format is the same 
as used in the program. 

1 . Store <x 0 , Co, and i Q . 

2. Store truncating values e and 

3. Store limit m to number of iterations. 

4. Store y. 

5. Leave space for a and ij. 

6. Set up storage for 10 Fourier series: X 4 , X 2 , 

X 3 , X 4 , W, v, h/h 0 , h n /h, n 0 Sz, and A. 

7. Store k 2 , k 3 , k iy A/2, and Bj2. 

8. Compute ij 2, sin (f 0 / 2), and cos (f 0 / 2). 

9. Store as follows: sin (-i 0 /2) — >X x 

0^X 2 
0^X 3 
cos (i 0 / 2)— >X 4 
O^TF 

0- *v 

1 — >ho/h 
1 — ho 
0 — >Uq5z 

0— s>A. 

10. Compute and store /3— f° 

1 + Vl-fio 2 

11. Compute and store the Fourier series 
l*=\ (l+^W) cos [F+( w )] 

+\ (1 — Vl — c 0 2 ) cos [F— (co)] — e 0 cos (w). 

12. Compute and store the Fourier series 
m*=\{ l+V^V) sin [F+(<*)] 

(1— Vl— Co 2 ) Sin [F— (w)]— e 0 sin (to). 

13. Compute and store 

=l-r+r 2 -r 3 + . . . +(— 1) V. 71=1, 2, 3, . . . 

l + r 

14. Compute and store the Fourier series 

f = : vdb? cos F +^ 2 cos 2F + ■ • • 

+j3 n cos 7 if\ 7i=l, 2, 3, . . . 


15. Compute the Fourier series product 

\p*=2 = [(XiX 4 — X 2 X 3 ) (X2X4T XiX 3 ) l*\. 

P 

16. Compute i/'* 2 . 

17. Compute 

a 0 _a 0 ( 1 \ / a 0 \ 2 ( a 0 Y 
P p Ki+ V )’ Vp/’ \p)' 

18. Compute the Fourier series (1— 3i p* 2 ). 

19. Compute the Fourier series (3 — 5 \p* 2 )\k*. 

20. Compute the Fourier series (35i^* 2 — 30)^* 2 +3. 

21. Compute 

h /aoV 

a 0 3 \ P / 

22. Compute 

/«o\ 3 (1 — 3^* 2 ). 

«o\ P / 

23. Compute 

*3 /aoV. 

«o* V P / 

24. Compute 

£=% Y (3— 5i^* 2 )^*. 
ffo \P/ 

25. Compute 

/ OnY. 

«o 5 \P/ 

26. Compute 

C=~ Y [ (35^* 2 — 30)\p * 2 + 3 ] . 

«o 5 V P / 

27. Compute S2*=A+i?+C. 

28. Compute -^p- by differentiating the Fourier 

series 0*. 

29. Compute 

Ao* 

pFA=_3A-4£-5C. 
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As the particle moves from 0 to P, r goes to r-t-dr. 
If we let i be a unit vector along r, and i-f-Ai 
be a unit vector along r + Ar we can draw the fol- 
lowing isosceles triangle, since both sides are 
unit vectors: 



We know from the definition of plane trigonom- 
etry, 

a- o ■ Ad 

Ai=2 sin — • 

Zi 

So we can write 


A i 
Ad 


2 sin 


A0 


Ad 


and since 


di Ai , 
d0 _ I™ A0 -1 ’ 


we see the magnitude of the vector di/dd= 1. 
Also clearly, from the diagram, as A0-O, Ai 
becomes perpendicular to i; so in the limit, 

lim unit vector perpendicular to i, 

dd a9->o A0 


which we will designate as j, and write 



Now, the radius vector r=ri, and the velocity vec- 
tor is simply 


But since 


d , . d , .. dr . . di 

TtV-mW-Tt'+'Tt 



di=jdd, 


and 


di_dd . 
dt~dt J ’ 


Then the acceleration vector a is dv/dt, or 

d 2 r . dr di (dr dd d 2 0\ . . dd dj 
& dt 2 1+ dt dt + \dt dt +r dt 2 ) i + r dt di 

However, using exactly the same reasoning as 
above, we can show that 



dj dd . 

dt dt *' 


So, substituting this expression for dj/df in the a 
equation and replacing di/dt with (dd/dt) j, then 
rearranging, we have: 

d 2 r dr dd .. (dr dd . d 2 6\ . dd dd . 

& ~df 2l+ didt i+ \dtdi +r dt 2 ) 3 r dtdt 
or 

a=[r— r(0) 2 ]i+[2r0+r0]j. 

However, 



so we can write 



where i and j are unit vectors along the radius 
and perpendicular to the radius vector, respec- 
tively. Therefore the component of the accelera- 
tion along r is 


f—rd 2 , 


and the component of the acceleration perpen- 
dicular to the radius vector is 


we can write the velocity vector 


v 


dr . . di dr . dd . 
dt' +r dt = di' +r dt > 


1 (L 

r dt 


(r 2 0). 
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30. Compute 

dO* 
d** = 


f-Y (3 — 15^* 
\a 0 / \p J do Vp / 


+ ^(^-Y (-60^*-fl40^* 3 ). 
«o \P/ 


31. Find 


dO dO* 


where the bar operation means that all F’s 
in the argument of the Fourier series 
become E’s. 


32. Find 

33. Find 


dO_dO* 
d^ — d^*' 

dO dO* 

d E~ d F' 

34. Find l— l* and m=m*. 

35. Compute M' , 


M' — — l 


+[©’ +,, (iu)] 2 ' ,co8E 

-(f) 


cos 2 E- 


1+r 


[e 0 2 cos (F+F)] 


+ 


[(2-Ol^+2(0]cos( F - E ) 
-(0 k cos (F-2E)] 

-[©’ + iu] e « cosf '' 

36. Compute N', 

"'“[©’-lU+fl 2to sin MrJ “ si " F 


— y sin 2 e +y^t v W sin (F+F)] 


(2— e 0 2 ) sin ( F—E ) 


1 + p 


+QJ [c 0 sin (. F-2F )]. 


37. Take the IF Fourier series stored, and com- 

, dTF 
pute w - 

38. Compute 

S=( l-e cos F) ^-(i+TF+^) e 0 sin f]- 

39. Compute 

X= (r) ( 1 +") 2 T *+ ~ cos E~[ 

\h 0 J L Vl-eo 2 /l-e 0 2 J 

40. Compute 

s '-(l57) A |>'( r 5 ?)+ m 'S] 


41. Compute 


S 2 =-=L= s. 

Vl-e 0 2 


42. Find the coefficients a, and a 2 of the sin F 

terms in <Si and S 2 , respectively. 

43. Solve for y in the equation 

a 1 +a 2 i/=0; (z/=-|) 

44. Construct 

dW 

m -Si+yS 2 . 

45. Compute 

cos %= — Xj 2 — X 2 2 +X 3 2 -t-X 4 “. 

46. Compute 

rri { d(j \ > / A \ dO 

1 =( — , ) COS ( r -A. 

\2Vl — Co 2 / \A>/ d^ 

47. Compute 

F 2 = T(— \ 3 l— X 4 m). 

48. Compute 

T 3 = T(X 2 l-\-\im). 

49. Find the constant terms in the Fourier series 

for Xi, T 2 , T 3 , and X 4 and call them b u b 2 , b 3 , 
and 6 4 , respectively. 

50. Compute 

b 2 i b 3 

a=~r and n= — y" 

01 0 4 
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51. Compute 

^=,=aX2+F(X 4 1 — X 3 m), 

dX 2 * | rp 

dE~~ aXl+T2 ’ 

^X 3 n i rp 

dE -^i+U, 

— v'^3 J rT(— Xj7+X 2 wi). 


57. Compute 

X>=( sin 2 + ^ + f) + ^ 1- 

58. Compute 

X 4 =(cos J+J~ f)+^. 

59. Store Xi and X 4 found in steps 57 and 58, along 

with X 2 and X 3 found in step 52, as the new 
X’s to be used in the next complete iteration. 


52. Compute 



where the tilde indicates that the constant 
of integration has not yet been added. 
We have the equivalence: 


X, — 5Xi, 


X 2 — ^X 2 , 


X 3 — 5 X 3 , 
X 4 =5X 4 . 


53. Compute C 7 =(Xi 4 -X 4 ) 2 d-(X 2 —X 3 ) 2 . 

54. Compute F=(Xi — X 4 ) 2 -|- (X 2 +X 3 ) 2 . 

55 . Find II and G, which are the constant parts 

of U and V, respectively. 


A B 

56. Compute and -=■ by iteration, where 
Z Z 


-H_( AX 

A 4 \2/ 

2 fo 1 * 4o 

cos 2 “+sm - 


and 


B 

2 



. IQ 

cos 2 — sm 2 


~ CdW 

30. Compute W— I (constant 


i n nlucl Arl \ 


of integration 


61. Operate to get W. 

62. Compute 

— U=i-w+(w)Y-(wY+ . . . . 

1 +W V / \ / 


63. Compute 


dn 0 8z 

~dE~ = 


1 +W 


^fF+Y)(l — c 0 COS E) 

(1 — r 2 )(l — e 0 cos E) 2 


y 


V 1 — 


e 0 


64. Find the coefficients a and /3 of the cos ( 0 ) and 

cos (E) 

65. Compute 


cos ( E ) terms, respectively, in 


A<7„=- 


-2a— eS 

2-e 0 2 ’ 


a n — 2/3-2c 0 q: 

A6l “ 2-r„ 2 

66 . Compute W= TF+ (ACo+ACi cos F) — store as 
new W. 


67. Compare a and (3 with e, and if |a| >e or |/?|>e, 
return to step 61 and use the W computed in 
step 66 for the W seen in step 61. 


68 . Compute n 0 Sz= 


J 


dn 0 8z 

dE 


integration is added. 


where no constant 


of 


69. Compute S, by splitting W series to get X, 
that part of W which contains terms inde- 
pendent of F, that is, select and collect all 
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Fourier series terms whose coefficients of F 
in the argument are zero. 

70. Compute Y, by writing TF=E+T cos F 

sin F, then letting F= 0, and subtracting 
X contained in step 69. 

71. Form S+eoT. 

72. Compute A by iteration: 

A=-| (E+e„T)+| (A 2 — A 3 + A 4 — A” + . . .). 

73. Compute (A+E+e 0 T)— store as new 

ho 

V 


74. Compute ^- = 1 + A — store as new 

75. Compute W by replacing F’s in step 66 by E’s. 

76. Compute v by iteration, where 

v=A— ^ (A+H- 7 )(l+r)— store as new v. 

A 

77. Compare 

\y n +i-y n \ with e, 

|«n+i — «»| with e, 

\Vn+i—v\ with e, 

and if any is greater, return to step 13 for 
next iteration. 
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